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centered-cubic crystal v/as developed to investigate the 
positional dependence of tlie diffusion coefficient. The 
simulation is an improvement over similar types of 
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simulation to investigate a more sophisticated system of 
diffusion, tlie face-centered - cubic crystal. The simulation 
provides o u t u t that can be used to d e v e 1 o n the diffusion 
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determine a functional form; for the diffusion coefficient if 
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INTRODUCTION 



A. !1IST0RY 

Since man first poured a colored liquid into water and 
observed that the liquid spread throughout the container, 
since he placed two metals together in heat and later 
observed that they melded together, or since he first 
noticed that the smoke from his fire mixed into the air, 
diffusion has been observed and its cause a source of 
wonder . 

Initial qualitative experiments in diffusion were 
conducted l)y Parrot in 1815 (Crank, 1975). In this 
experiment, he noted that no matter hov/ carefully he 
controlled mechanical agitation and convection, the separate 
gases in the experiment always diffused throughout the 
container. One other important precursor development was 
Fourier’s development of the equation for heat conduction 
(Fourier, 1822). 

The history of the scientific study of diffusion is 
rooted in a paper by Adolf Fick (1855). Fick used Fourier's 
mathematical methods to describe the equations of 
diffusion. One equation states that the rate of diffusion 
is proportional to the concentration gradient across the 
[)lane of diffusion: 

J =- DO C/ 8x), (1) 
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where J is the flux of diffusing material, D the (i i f f ’i s i o r. 
coefficient, and C the -concentration of the diffusing 
naterial. The units of the diffusion coefficient are 

(length)'^/(time), and it is usually given in centineters 
per second. Equation (1) is known as Pick’s first law of 
diffusion. In this v/ork, Equation (1), is assumed to be of 
this form, even for non-constant D. 

Pick also developed a second equation \^fhich describes 
the time rate of change of the concentration of the 
diffusing material. Known as Pick’s second law, it states 
that the partial derivative of the concentration v/ith 
respect to time is equal to tlie negative of the divergence 
o f the f lux : 

9C/ at = - V J . ( 2 ) 

In one dimensional cartesian coordinates, this reduces to: 

aC/ at = - 3J/ 3x . (3) 

\vhen Equation (3) and Pick’s first law (Equation 1) are 
combined Equation (3) becomes: 

ac / 3t = { a / ax [ ( D ) ( ac / 3x ) ] ) , ( 3a ) 

and, if one assumes that the diffusion coefficient is 
constant with respect to posit i. on, the equation further 
reduces to: 

ac/a t = D( a^c/ ax^) . (3b) 



P. 



The general solution to this equation, given an initial 
impulsive planar source h i c h then spreads across a o n e - 
dimensional infinite material, is: 

C=(A/t^^^)exp(-x~/4Dt), (4) 

where A is a constant, t the elapsed time, and x the 
position of interest within the material. A consequence of 
this solution is'that after sufficient tii"e Liie 
concentration will become uniform throughout the material as 
seen in Figure 1. 



C 




Figure 1. Concentration Profiles at Different Times 

Einstein (1905) provided an atomistic definition of the 
diffusion coefficient: 

D = <j<^>/2t, 

where 21 is the vector displacement of a given atom from its 
initial position to its final position at time t. The Dirac 
1) r a c k e t. s denote an average over a large number of diffusing 
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atoms. It can readily be shown that t h e <1 e f i n i t i o n s for D 
in Equations (1) and (5) are equivalent. 



B. DIFFUSION MECHANISMS 

Since Pick's and Einstein's v/ork in diffusion theory, 
many people have tried to explain more fully tiie mechanisms 
by v/hich an atom diffuses through another material. I'r. e 
majority of this work was done by Kirkendall (1942), Hartiy 
(1946) and Darken (1948). These men developed tools to 
explain diffusion based on a knowledge of the nature of the 
atom and the method by which atoms combine into bull: 
materials. Given the concept of the crystal structure of 
solids, many mechanisms of diffusion can be identified. 
Table 1 is a partial, list of the mechanisms of diffusion 
through solid crystals: 



TABLE 1 

PARTIAL TABLE OF DIFFUSION MECHANISMS 

1. Interstitial mechanism 

2. Inter stitia Icy mechanism 

3. Crowdion mechanism 
4 . Vacancy m e c li a n i s m 
5. Relax ion mechanism. 



The INTERSTITIAL M ECnANIS M d e s c r i 1) e s 
small atoms through a crystal structure of 
atom advances through the structure by 



the 


diffusion 


0 f 


larg 


e r atoms. 


T h G 


m 0 V i 


n g from 


0 n e 
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vacant L n t e r t i u i a 1 site 'a' i t h i n the crystal to a n o i h e r . This 
is s h o v; n in Figure 2 . 



Figure 2. Intertitial Mechanism of Diffusion. 

V.’hen the interstitial atom's size approaches that of the 
lattice atoms, they move by the INTERSTITIALCY . 1 E C II A N ISM . I n 
this method the defect atom moves into a normal lattice 
position by pushing the lattice atom out into an 
interstitial site. Figure 3 illustrates this mechanism. 





Figure 3. Intertitialcy Mechanism of Diffusion. 

The C R 0 W DIO N' M E C 11 A N I S M of diffusion requires a defect 



atom with an energy large enough such that it can push its 




way onto a lattice site, compressing a close packed rov; of 
atoms. This mechanism can be thought of as a line of 
billiard balls touching each other. Then a cue ball liits 
tills line so that each ball moves one position, replacing 
its neighbor. An example is shown in Figure 4. 




Figure 4. Crowd ion Mechanism of Diffusion. 



All crystals have vacancies in their lattices at a 
finite temperature. In the V A C A N C Y MFXHANISM of diffusion, 
these vacancies are filled by the diffusing atoms. This 
method of diffusion is predominant in self diffusion of a 
material since a lattice atom and the vacancy simply 
exchange position. Figure 5 illustrates the method by v/hich 
the vacancy moves. 




Figure 5. Vacancy Mechanism of Diffusion. 
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The R E L A X 1 0 X M E C H A X I S !I of diffusion is a n o d 1 f i e o 



x’acancy mechani-sn. In the area of a vacancy thc-> lattice 
atoms relax tov/ards the vacancy. These atoms, and any defect 
atoms in this area, are able to move by an irregular 
movement similar to motion of a liquid. This mechanism can 
be seen in Figure 6. 




Figure 6. Relaxion Mechanism of Diffusion. 

J.R. Manning (1968) has an excellent discussion of the 
mechanisms of diffusion. 



C. BASIC CONCEPTS 

Crystaline solids are made up of atoms that arranp, e 
themselves in various lattice structures. There are seven 

these are: 



3 



types of crystals. 



TABLE 



/ 



TA]]LE OF CRYSTAL STRUCTURE 



1 . cubic 



2. tetragonal 



3. orthorhombic 



4. hexagonal 



5. rhombohedral 



6. nionoclinic 



7 . triclinic 



I'/ithin each of these types of crystals subdivisions exist 
that further divide these seven structures into a total of 
fourteen crystal structures. L.A. Girifalco (1964) gives a 
good discussion of the different crystal structures. The 
face-centered cubic (fee) crystal structure is used in this 
thesis . 



Figure 7. Different Propectives of the fee Crystal. 

Figure 7 illustrates different sections of a fee crystal. 
There exists sites between the atoms of the structure known 
as interstitial sites, see Figure 8. 





(1 0 0 ) 



(3-D) 
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0 -interstitial sites 



^ -lattice atons 



Figure 3. Tntertitial and Lattice Sites of fee Crystal. 

Each site in the fee structure, whether an interstitia] site 
or lattice site, has t\\felve nearest neighbors. Figure 9 
shows the relationship of a site and its nearest neighbors 
within the crystal. The interstitial method of diffusion 
described above moves an atom to one of these nearest- 
neighbor sites. 




0 - interstitial sites 




lattice sites 



Figure 9 . Fee Atom and its Twelve e a r e s t Neighbors. 
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The force acting between a pair of atoms changes as the 
distance i)ctu'een the atoms changes. Initially, at largo 
separations, the force is attractive, becoming stronger as 
they approach one another. At small separations the force 
become repulsive. This force is best described by a 
potential energy function. The minimurn of this potential 
function is the equilibrium position between the pair of 
atoms. Figure 10 shows a representative plot of tliis 
potential energy for a pair of atoms. 



Repulsive 

Attractive 




Figure 10. Representative Potential Energy Plot. 

In a crystal structure this potential energy function is a 
multi-body problem. The result of tliese extra bodies is 
tl'iat at a greater distance, the function oscillates about 
the zero energy position as shown in Figure 11. 
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Figure 11. Energy PoEentia] of a Hulti-Body Systen. 

An atom has a kinetic energy that is characteristic of 
the thermal energy of its surroundings. This motion is 
known as the thermal motion of the atom. 

When the thermal motion and the interatomic potent:: a 1 
of the crystal structure are combined, the atom's motion 
becomes a vibration. An atom oscillates about its site, 
always moving near its lattice position. Add a diffusing 
atom, in an interstitial site of the crystal, witi: its o\vn 
vibration and its own potential energy and this interstitial- 
atom will gain and lose energy as it vibrates between the 
lattice atoms. Figure 12 shows a qualitative graph of 
energy versus time of an atom's history while in an 
interstitial site. 
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The potential energy of the atoms in a crystal structure 
sets up an energy barrier that a diffusing atom must 
overcome before it can move to one of its neighboring sites. 
This energy is known as the activation energy for diffusion. 
Figure 13 shows a qualitative graph of an atom’s kinetic 
energy with respect to time, with the activation energy for 
a move indicated. 



KE 




Figurel3. History of Atom’s Energy Noting Atom's 
Activation Energy. 
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D. ra::dom walks 

Early in the study of diffusion n e c h a n i s m s , Che 
possibility of simulating the overall motion of diffusi. on 
was investigated. The method devised is still used today 
and is known as a Random V/alk Process. A random walk in one 
dimension can be described as assigning a direction, either 
forward or backv/ard, to the faces of a coin. \vhen the coin 
is flipped, one step* is made backward or forward according 
to the face of the coin displayed. Let be the vector 

connecting the origin to the final position of the atom. 
The equation of R^ is then: 

R^ = r^ + t2+...= 

where rj^ arc the steps representing the indi. vidual jumps of 
the atom. If each direction of motion is equally likely anci 
individuai jumps are of equal magnitude, then the magnitude 
of R^ is 

Low let a number of atoms, initially at x = 0, each complete 
one dimensional random v/alks. Close to 95% of the atoms 
should be within + or - of the origin. 

E. CCdlPUTER SinULATION OF DIFFUSION 

Shortly after the development of the electronic computer 
liayer and Ulam proposed that computers be used to simulate 
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physical sysf. ems, !. etropolis (1953), and Kin 3 (1951) 
suggested that the lionte Carlo method of simulation (this 
method v/ill be discussed in depth below) could be used to 
simulate the random walks of particles. The problem was 
that, at the time of King's paper, the computational 
capabilities of the machines in use were easily exceeded by 
Monte Carlo simulations. By 1951, computer technolog 3 ' was 
sufficiently advanced that Flinn and McManus were able to 
complete a Monte Carlo simulation of solid state diffusion. 
However, as late as 1970, there still had been no simulation 
of the diffusion processes based on the random walk theory. 
By 1970, it was evident that a new simulation process was 
required since the analytical processes of diffusion were 
unable to solve the topical diffusion problems evolving in 
the literature. In 1971 Bennett and Alder published the 
results of their simulation of diffusion by vacancies and 
divacancy random v/alks. Their paper marked the beginning of 
the random walk simulations of diffusion. Through, thc: 
1970's and 1980's, the field of Monte Carlo simulation of 
the diffusion process came into its own. Simulations have 
been used to test the correlation factors of the diffusi. on 
coefficient (Guy et al, 1977; Guy, 1973), the self diffusion 
of tracer atoms in an alloy (B.akker, 1976), and many other 
specific diffusion mechanisms (Murch, 1984). 
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CONCENTRATION DEPENDENT DIFFUSION COEFFICIENT 



Concentration dependent diffusion v;as' apparent oarl}' in 
the study of diffusion and the diffusion coefficient. C. 
M a n t a n o ( 1 9 3 3) made the first attempt to m a t ’n e m a t i c a 1 1 y 
solve for the experimental results that showed the 
concentration dependence. Using the experimental results of 
metal - metal diffusion, lie found that the concentration 
dependent diffusion 'coefficient would have the follov/ing 
f or m : 



D(C) = (-l/2t) (dx/dC) / xdC. (8) 

This equation is known as the Matano solution for the 
diffusion coefficient. 

Bowker and King (1978) attempted to develop a Monte 
Carlo simulation to determine the concentration dependent 
diffusion coefficient. Their first attempt used a two 
dimensional lattice gas model of a simple square plane and 
their results appeared to correspond with the ex peri, mental 
data of oxygen-tungston diffusion. However, a more 
detailed investigation of the simulation led to discovery of 
problems in the data collect i. on of the concentration (Murcli, 
1981). M u r c h (1980) also tried to develop a simulation Lliat. 
could provide a concentration dependent diffusion 
coefficient (this simulation is discussed belov/), but his 
simulation showed no concentration dependence. 



A recent paper by A . J . S 1 a v i n and ? . I' . U n d e r ii i 1 1 
( 1 9 A. 4 ) a 1 1 e n: p t e d to prove t b. a t a concentration d e p e n d c n t 
diffusion coefficient could not be found using a randon v/alk 
simulation of diffusion. However, a paper by H. Ghez and 
vJ.E. Langlois (1986) showed that a concentration dependent 
simulation can be developed if care is taken in the 
collection of the fluxes at the midpoint betv/een the planes 
of the crystal. 
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OBJECTIVE 



II . 

The {) r i HI a r y objective of this thesis is to i n v o s t i a n t e 
the possibility that the diffusion coefficient is position 
dependent. The idea that such a dependence night exist is 
not new. In fact, a positional dependence has been observed 
in a crystal v;hen there is a tenperature gradient across the 
structure. It has always been accepted t 'n a t w h'e n 1 1 * i s 
gradient is absent, and the temperature is const, ant 
throughout the structure, the diffusion coefficient would be 
independent of tiie position of the diffusing atom. The 
idea of a position dependent diffusion coefficient \;as 
proposed by Collins (1985) as one way to test the assumption 
that the partition function, in Gibbsian statistical 
mechanics, does not divide into purely positional and 
velocity parts. 

Though proposed early in the history of simulation, 
diffusion simulations did not actually begin until Bennett 
and Alder (1971) developed t li e i r program. They were studying 
persistence effects in hard sphere gases and used a simple 
ilonte Carlo simulation to find the vacancy correlatio:i 
factor in a two-dimensional fee crystal. The first actual 
simulation of three di. mensional diffusion was done by 
d e B r u i n and 1 i u r c h ( 1973). This simulation attempted l. o 
develop a method to determine the diffusion correlation 



factor. 7 lie diffusion correlation factor (f) is based on 
the Einstein diffusion coefficient method of determining the 
coefficient. Redefining Equation 5, 

D=nX^f/2t (9) 

the diffusion corellation factor can be defined as 

f ■= <x2>/nx2 (10) 

where n is the number of jumps the diffusing particle tal:es, 
X is the distance of each jump, and <x > is the mean-square 
displacement of the diffusing particle after time t (Le 
Claire, 1970). Early simulations of diffusion dealt with 
either diffusion of a vacancy, or self diffusion, and all 
were based on the Einstein diffusion equation (see Equation 

5). 

The first reported use of Pick's First Law in a 
diffusion simulation was by Guy et al. (1977) and Guy (197S) 
when the correlation factor was calculated. The method used 
in these investigations was to place an initial 
concentration of vacancies in a chemically homogeneous 
square planar structure and then allow the vacancies to 
move throughout the structure by a random v/alk process for a 
period of time. This method was found to have problems in 

that they had assumed that the atomic flux and the vacancy 
flux across the plane were different, and set up the 
simulation in that vein. In 1979 ‘d u r c h and Thorn again 



24 



tried to simulate the tracer correlation factor by this 
method-. Their method v/ a s to set up a t h. r e c - d i m. e n s i o n a i 
array with an inexhaustible reservoir of tracers or 
vacancies at one end (x=0), and at the other end (x=l) 
remove the tracers into a sink. The array was periodic lii 
the y and z directions. From this simulation one could 
calculate the net fl. ux across a plane, and from this 
calculate the correla'tion factor. In their paper, iiurch and 
Thorn reported that this method gave as accurate answers as 
the Einstein method and was not so difficult to use. 

Murch (1980) also used the Einstein method (Equation 
5) to test the effect of diffusing atoms concentration on 
the chemical diffusion coefficient. This simul.ation was 
based on a three-dimensional simple cubic array. Within this 
array, he placed two reservoirs of atoms eight lattice units 
in from the x boundaries, each with a separate chemical 
potential. He then allowed the atoms to diffuse using a 
nearest-neighbor interaction, interstitial mechanism. lie 
concluded from those investigations that the concentratic^n 
did not affect the diffusion coefficient. 

Pick's second lav/ was not exploited in a diffusion 
simulation until 1978, v;hen Howker and King used it to 
study the diffusion coefiicent. After initializing a lattice' 
of vacancies, they placed a concentration of defects in one 
half of tlie lattice, and a different concentration in the 
other half. After a number of random steps of those 
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defects, the concentration across the structure was counted 
plane by [ilane and a concentrati. on profile was flev elope cl. 
Froia this profile the diffusion coefficient was derived 
n u n e r i c a 1 1 y . M u r c h ( 1 9 <3 1 ) used the second lav/ n e t h o d but 
with the external reservoirs that he and Thorn (1979) 
developed for their first law investigation. He noted that 
if he used different data planes to develop his 
concentration profile’ he found different results. This led 
him to believe that the method had some problems with the 
development of the concentration gradient across the 
structure. 

The development of these simulation techniques lias 
proven that the diffusion coefficient can -be studied using a 
digital computer. Each of these simulations used the 
premise that the coefficient is constant throughout the 
structure, but M urch (1980) suggested that these techniques 
could be used to see whether the diffusion coefficient has a 
positional dependence. Another limitation to all of tliese 
simulations is that the authors used either the Einstein 
method to determine the diffusion coefficient or Pick's 
first or second law. The simulations were never analya. ed 
v/i. th more than one technique. 

A comparison of the diffusion coefficient as calculated 
by Pick's first and second laws will be the primary ni c t h. o d 
to determine the positional dependence in this thesis. If a 
positional dependence is recorded, a method other than 



ride's laws will be necessary to determine tlie diffusion 
coefficient. This method i, s !c n o w n as the ' F o k k c r - P 1 a n !•: ' 
equation of diffusion: 

8c /9t = -8 /8x[A(x,t)C(x,t)]+ 1/28^ /8x^[D(x,t)C(x,t)] (11) 

In this Equation (11) the first clement of the sun 
constitutes a drift factor where A(x,t) is a drift velocity, 
and the second element is the diffusion equation. The 
equation actually is giving the time evolution of the 
probability function, C(x,t). In the simulation this 
probability function is the concentration. If the time 
derivitive, spatial derivatives, and the functional form of 
the diffussion coefficient are known, the diffusion 



coefficient can be fit to the simulation data. 



III. COMPUTER :;ODEL AMD SIMULATION 



A. PHYSICAL MODEL 

A computer simulation model must mimic a physical model 
or experiment. In this simulation the physical model, is a 
simple face-centered-cubic crystal membrane which is 
infinite in the x and y axis and very thin (less than 50 
lattice planes) along tlie z axis. The mechanism \;ill be 
interstitial diffusion. This means that the diffusing atoms 
are very small in comparison to the size of the lattice 
atoms. These diffusing atoms start in an inexhaustible 
reservoir outside the input plane of the crystal, and are 
completely absent from the volume outside the other'side of 
the crystal at time equal to zero. At this starting time, 
the atoms begin to diffuse. After a period of time the 
distribution of diffusing atoms approach a steady state; in 
other words the transient effects on the concentration 
profile have ceased. At this time the flux across onc!"! plane 
in the crystal structure is checked and the concentration 
gradient is recorded. The diffusion coefficient is then 
found from these values. This thesis c o n s i. s t s of a c o Ci p u t e r 
simulation of this physical s y s t e tn , v; t. i c ‘n was designed to 
study the diffusion coefficient. 
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7 ;i e p 'n y s i c n 1 r;. o d g 1 is t r n :i s ,i g t e d in s u c h a a y 1 1 1 r. t t ii e 

c 0 m j; u t e r can f o 1 1 o the actions t li a t take place. ( i i s j s 
I: n o n as c o ;n p u t e r s i n u 1 a t i o n or c o r. p u t e r p. o d o 1 i ;; o i a 
p h. y s i c a 1 s y s t c p . ) C o n: p u t c r simulations can be separatee' 

into tv;o families: continuous-time simulation and c'iscreto- 

event simulation. Continuous-time simulati. ons (also !:no\;r. 
as t i :n e s t e p simulations) simulate events w Ii i c ii can 1. e 
described by s i m u ]. t a n e o u s equations. The e (i u a t i o a. s 
describing events arc a], lov/eei to proceed for a f)re- 
detorminedperiod of time (a time step), and tlien infor::iation 
is placed back into the equations which are to be s o I e d. 
numeric a ]. ly. -Discrete-event simulations are used for 
pliysical systems w'hich cannot be described by a m a t ’nc m. a t i ca 1 
function. The simulation runs from= event to event ra. thor 
then advancing a number of actions during a time step. 

Simulations can also be divided into typos by t a. e 
!ii c t h o d s used to make decisions. If the simulation uses 
r a n (! o !ii numbers to m a I: e decisions, it is called a i . o n t c C a. r 1 o 
simulation. This is because early simulations (those prior 
to computers) used devices, such, as dice, to c ’n o o s e t hi e 
ran do. m numbers that were used to ma!ce t!ie decisions. Once 
simulations became more sophisticated, -the community tried 
(in vain) to shift to the name 'stochastic simulation', but 
the name ' i o n t c Carlo' has stayed popular. In t ’n e o t ii o r type 
of simulation, called 'deterministic s L ::: u 1 a t i o n ' , a 1 1 



decisions are made by rules set down prior to the actual 
running of the sinulation program (Harr is on, 19 8 5). The 
simulation used in this thesis is of the ’lonte Carlo type. 

C. RANDOli WALK 

The diffusion events are simulated by a method called a 
three dimensional random walk process. This process is 
simple in form and gan easily be set up. Consider a two 
dimensional example: Imagine that a marble is on a Chinese 

Ciieckers board. The marble has six neighbor sites to v/b. ich 
it can be moved. A random \'/alk of the marbles can be devised 
where each direction that the marble can move is assigned a 
number between one and six. A die is rolled and the marble 
is moved in the direction indicated by the number on t'ne 
face of the die. This is an example of an uncorrelated two- 
dimensional random walk. Other forms of random walks are 
the 'v/eighted’ and 'correlated' random walk. Returning to 
the Chinese Checlcer board example, if you move tv;o marbles, 
one that is black and one that is white, with the sam. e roll 
of the die, that is, if every time the die is rolled for tb. e 
white in a r b 1 e you a 1 s m. o v e the black, marble, then the t o 
moves are connected atid the random v/alk is correlated. ;\n 
example of a v/eighted random walk can be seen by t a !c i n g t \/ o 
dice and assigning one number to each of the four 
directions you want to have a lesser chance of the narl)1e 
moving to, and assigning three numbers to each of t'ne other 
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t v; o directions. T ii i s is a e i g h t e d r a n d o n w a 1 '< in the 
direction of interest. Tiie random v/alk used in this thesis 
is a n o n - w e i 3 1 1 1 e d , u n c o r r e 1 a t e d random walk. 

D. COMPUTER CODE 

The code for the computer simulation is divided into two 
separate programs, ^ DIFSET ' (for diffusion simulation set-up 
program) and ' DIFFUSE ' (for diffusion simulation program). 
These programs combine to form, a Monte Carlo simulation of 
an interstitial mechanism type of diffusion. 

1 . Dif set 

The program 'difset' can be separated into three 
[)arts, each of vvMnich intializes a table: 

TABLE 3 

TABLE OF DIFSET ARRAYS 
Crystal site table 
Nearest-neighbor table 
Site vs plane table. 

a. Crystal Set Up 

The face-centered-cubic crystal used in tiis 
simulation is set u[) in a three dimensional array called 
CRYST3(I,J,K). The three dimensional array variables must 
meet certain limitations. T!i.ese limitations are that the T 
and J variables have to be even numbers. There are no 
limitations on the K variable. Tlio reason for thcise 
limitations is that the T and J axis in the simulated fee 
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crystal are set up to have periodic boundaries. The 
initial array is a perfect face-centered-cubic crystal v. ith 
no defects or vacancies. An array position designjites the 
center of a site and gives no indication of the size C)f an 
atom that might occupy the site. Atom size is not a factor 
in the simulation or in the random walk used, and therefore 
is not included in the set-up beyond the' assumption that the 
diffusing atoms are 'small enough (relative to the lattice 
atoms in the crystal) to allow the interstitial mechanism of 
diffusion to occur. The set-up of the array CRYST3(I,J,T) 
is accomplished by placing the number tv/o (2) in the array's 
register as an indication that a lattice atom is in t'aat 
array position. To accomplish this, all the array positions 
v;here the I, J a n<l K variables add to an odd number are 
assigned the number two (2). Table 4 lists examples of the 
initial value set up in 'DIFSET' for the lattice positions. 



TABLE 4 



EXAMPLES OF ARRAY ENTRIES TO THE CRYST3 ARRAY. 

CRYST3(3, 1 , 1) = 2 
CRYST3C 2 , 2 , 1 ) = 2 
CRYST3( 2 , 1 , 2) = 2 
CRYST3(5 , 2 , 4) = 2 



The interstitial 



positions of the crystal 



are set up in the 



same way. 



Site 



indices that add to an 



even number are 



assigned a zero value to indicate that the interstitial 
p o s t i t i o n s are e t'. p t y . T 'n e crystal's first plane 
interstitial sites are filled by diffusing ato.ns. Tl'ierefore 
the array positions (even,odd,l) and (odd,even,l) are given 
the value of one (1) to indicate that these interstit. ial 
sites are filled with tlie diffusing atoms. The array is no\-/ 
set and an example plane is sho\i/n in Table 5. 

TABLE 5 

EXAliPLL OF THE FIRST TUO PLANES IN A 4X4XK CRYSTAL. 

CRYST3(I , J , K) plane //I 

2 12 1 
12 12 
2 12 1 
12 12 

CRYST3(I, J,!C) plane #2 

0 2 0 2 
2 0 2 0 
0 2 0 2 
2 0 2 0 



The three dimensional array that was set u[) is 



difficult for the 


computer 


t o 


use, so a 


one dimensi 


0 n a 1 


array, 


C R Y S T ( L ) , 


is set 


u p 


to be e (| u 


ivalent to 


t h e 


CRYST3(T 


, J , K ) three 


dimensional 


array. For 


methodology 


see 
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code i-n Appendix j3 . 



b . ! ! e a r e s t - ! I e i g ii b o r Table 

T he next n a j o r step of the ' D I F S E T ' p r o g r a r. is 

to set up a nearest-neighbor table that \;ill be used to 

simplify the random v; a 1 ic process. To r e v i e v; t !i e diffusion 
by interstitial method, there are twelve interstitial sites 
v/hich are nearest-neighbors to any other interstitial site. 
Therefore, a diffusing atom must move into one of these 
t V/ e 1 V e n e a r e s t - n e i g li b o r sites. The array C R R K B ?. ( L , E ) (for 
crystal nearest-neighbor array) is set up such that the 
value held in the register is the L variable of anotlier 
interstitial site in the CRYST(L) array. Thus if you \'/ant 
to knovv^ the twelve nearest-neighbors to CRYST(13), you 
would look them up in the CRMRBn(13,R) where the t \/ e 1 v e R 
values are the twelve nearest-neighbor site L variables. 
The periodic boundary conditions for the I and J directions 

are set up in this array. The boundaries in these 

directions are made periodic by setting up the neighbor 
table to take an atom that leaves the boundary of t ii e 
crystal on the left side, and replace it on the right 



a e 



of the 


crystal as if the 


left and ri 


gh t 


side were 


neighbor 


planes 


. The edge sites of 


t ii e K a x i s 


are 


the entry 


and exit 


p o s i t i 


0 n s of diffusing 


atoms and 


a s 


s u c ii have 


special 



n o a r o s t - n e i g h 1) o r table values. The entry edge is simulated 
such that the plane is i. n contact with a infinite supply of 
the diffusing atoms, so the interstitial sites arc always 
filled with the d i f f u s I n <’ atoms. A n o t Ii e r way of looking at 



t h i sa s s u iri p't i o n physically is: if a diffusing a t o leaves a 



surface i n t e r s c i L i a 1 position, another d i f f u s i. n g a L o ir. L a e s 
the interstitial position before the simulated times top is 
completed. T It ore fore, the neighbor table is sot. up to 
contain only the four neighbors that are v/ithin the crystal. 
This means that a diffusing atom has the opportunity to 
enter the crystal only one third of t.he time. An example of 
this table is located, in Appendix F. 

c. Site vs Plane Table 

The final item that the program 'DIFSET' 
develops is the lattice site plane table. This array 
ZPLM(L) (representing the plane ft in ythe z direction) 
correlates the site number L to the plane number v;iiich is 
contained in ZPLM(L). 

The information developed in 'DIFSET' is then 
placed into an output file to be used as input to the other 
program 'DIFFUSE'. 



2 . D 1 f fuse 

'DIFFUSE' is the actual simulation prograns that 
models the diffusion process. This is done using a Konte 
Carlo timestep simulation. At each tin estop, th. e 
simulation lool: s at each crystal site. If it contains a 



diffusing 


a t o n 


the 


simulation 


then 


1 ooks 


; a t a 


nc ip 


libor site 


selected 


at r a n d 0 


m. If 


the 


no ig 


h b o r 


site 


is e 


m p t y t li e 


diffusing 


atom 


i s 


moved 


into 


that 


no ig 


h b o r 


site 


. In 1 1. e 
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simulation, the process is cl o a e v/ i t h two tables of rancior'- 
numbers. At the beginning of each t i n e s t e p , a random n u n ’n e r 
generator is colled to make a table of random numbers, one 
for each site in the CRYST(L) array. The random number is 
uniform, 1 to LMAX, where Li'lAX is the total number of sites 
within the array. The sites are characterized in this way, 
so the’ simulation can be us'ed later to study diffusion by 
the vacancy mechanism with little modification. .\n other 
table of random numbers uniform 1 to 12 is set up. Again the 
size of the table is ecjual to the number of sites in the 
array CRYST(L). As a time saver, during the first few 
timesteps the number of random numbers generated is equal to 
the number of sites in one plane multiplied by the timestep, 
and only those first planes are used in tlie simulation. 

After the random numbers for the timestep are 
assigned to their tables, the simulation calls the first 
random number (from the uniform 1 to LMAX set), sets that 
number equal to L and looks to see if CnYST(L) contains a 
diffusing atom. If the site does not contain a diffusing 
atom, another random number is selected from the first 
table. If there is a diffusing atom in the site, a ranciom 
number is selected from the second table. This number is 
set equal to R and the nearest neighbor tab-le is used to 
find the value in CR liK RR( L , R ) • This value is then set equal 
to LI and if CRYST(Ll) equals zero tiie diffusing atom is 



r.ioved into the new location. This process continues until 
all L ! i A X numbers 1 1 a v e been tested. 

During each of these moves a number of factors 
required for the n u m e r i. c a 1 solution of t li e diffusion 
coefficient are calculated. 

a. Concentration 

The concentration of atoms in each plane is 
determined at the beginning of each timestep by looking at 
every site in the array. If the register contains a one ( a 
diffusing atom), then the concentration in that plane is 
increased by setting K = ZPLN(L), and adding one to t'ne 
value already in the array CONS(TilST?,k). The concentration 
is then divided by the number PLMMAX, v/hich is the total 
volume of a single plane. This puts the concentration into 
units of particles per unit volume. These units are needed 
for use in Pick’s law. This value is stored in. the array 
CVOL(TMSTP,K) . 



b . P 1 u X 

The flux density through each plane is 
determined when an atom is being moved from the old site to 
its new site. Since the flux density is actually bet\/een 
two planes, by definition the area bet\v'een plane K(l) on the 
left and K(2) on tlie riglit will be considered part of plane 
u(l). Th. is means that if an atom diffuses from K(l) into 
K ( 2 ) then the flux of K ( 1 ) has 1 / P L X M A X added to it. 
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II o V.' e V e r , if the atom moves from K ( 2 ) hack to II ( 1 ) the flux 



of II ( 1 ) is decreased by 1 / P L >' 1 1 A X . T h u s i L U X ( T 1 1 S T P , X ( -1 ) ) ) = 
F LU X ( T II S T? , II ( 1 ) ) + 1/PLXMAX if the atom moves from X(l) to 
11(2) and FLUX(fl!STP,K(l)) = FLUXCTIIST? , K( 1 )) -1/PLII' AX if it 
moves from II ( 2 ) to K ( 1 ) . 

c. First Derivative 

The first derivative of the concentration v/it!i 
respect to position, is calculated by subtracting the 
concentration per unit volume of the plane after the plane 
where the derivative is required, from the plane preceding 
the plane where the derivative is required. For example 



DERC1(TMSTP,II) = CVOL(TMSTP , K + 1 ) - C VOL ( T M STP , II - 1 ) . The 

units for the concentration gradient, as stated above, are 
particles per inter-planar distance to tiie fourth po’u’or. 
Ivith these units the change in tlie position is 1 and the 
mathematics is simplified. 

d. Second Derivative 

The second derivative v;ith respect to the 
position is obtained in the same manner, w i t li the exception 
that the 1st derivitve is used rathier than the planar 
concentration. D E R C 2 ( T li S T P , II ) = D E R C 1 ( T I ! S T P , II + 1 ) - 



D L R C 1 ( T I ! S T P , K - 1 ) . The units for the second derivative are- 
particles per in tor- planar distance to the fifth j'ov/er. 



are 



derived in 



e . Time Derivative 



The time derivative 



q r, -- 



manner as the positional derivative. Ti.e concentration of 
the plane at the beginning of the timestep is subtracted 
from the concentration at the end; BERT(T!iST?,iC) = 
CVCL(T:iSTP+l ,!•:) - CVOL(Tr:STP,K). The units for this factor 

are particles per tin estep per inter-planar distance cubed. 



f . Averaging 

After a pre-de ter mined number of times tops iiavc 
been completed, t'ne arrays containing the information are 
averaged. This gives another set of arrays that are used to 
derive ’the diffusion coefficient. The single times tea 
arrays are tlien reset to zero, and another set of tin estops 
are begun. 



g. Diffusion Coefficient 

In this simulation the diffusion coefficient is 
derived by two methods so t ii a t the results can be c o r: par o d 
and the [) o s i t i o n a 1 dependence of the coefficient may be ~i o r e 
easily seen. The first method uses F i c ic ’ s first 1 a . 
Having found the average flux over a number of tiinesteps, 
A FLUX ( ATi 1ST? , H ) , and the average first derivative. 



.'\ jF.RC 1 ( ATiiST P , X ) , the diffusion coefficient can be found 



b 



dividing the average flux by tlie negative of the derivative, 
D I F C 0 1 ( A T ;■! S T P , K ) = - A F L U X ( A T 1 ! S T P , K ) / A D F H C 1 ( A T i i S T P , K ) . T Ii e 
diffusion coefficient is not determined for everv timesten 



O (\ 



because, at t i ::i e s , t !i e first derivative for a s i- n 1 e 
t i p. e s t e p is zero and the r e s u 1 1 i ti ^ d i f f u s i o n c o e f £ i c i e 1 1 1 is 
infinite. This result v o u 1 d s k e v; t it e a v e r a e e d diffusion 
coef f icient . 



The second method of determining the diffusion 
coefficient is based on Ficic's second law of diffusion. 
More specifically, the second law holds that tiie 
diffusion coefficient- is constant with respect to position, 
9c/ 8 t = D8^c/ 9x“. T!ie diffusion coefficient is then 
determined by dividing the time derivative of thie 
concentration by the second positional derivative, 
D I F C 0 2 ( A T : ! S T P , K ) = A D ?• P/f ( A T : I S T P , k ) / A D F R C 2 ( A T ! ( S T P , k ) ) . 
Averages are used here for th. e same reason as described for 
the first method. 
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T’v . RESULTS 



A. PROOF OF SIMULATIO:: ' S FEYSICAL ECU I V VLE'XE 

The proof of the simulation's ability to model the 
physical v;orld can be seen in three v/ays. First the 
simulation adheres to the conservation of matter. Seconr’, 
tlie diffusion coefficient that is developed from the 
simulation's data lies within the limits of t'ne 
experimentally found values for • the diffusion coefficient. 
Third, the time average of tiie concentration profile 
approaches the infinite time result predicted for the 
concentration . 

1 . Conservation o f r i a 1 1 e r 

The conservation of natter and Pick's second la;: 
(Equation 3) are equivalent. Therefore, if Equation 3 lioltls 
in the simulation then the conservation of matter holds oy 
definition. To tost that this requirenont is met at all 
times during the simulation, the first derivative of t i i o 
flux density \/ a s determined and compared to t ii e t i. n e 
derivative of the concentration during one run of tiio 
simulation. It v/as found that tlie conservation of matter 
held for all time in the simulation. 

2 . Comparison o f the Diffusion Coefficient 

The diffusion coefficient developed from the 
simulation data is in units of inter-planar distance squared 
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p o r U i r: e F I o p . In order L o c o p a re t d c d i i' f u s i o n 
c o e i r i c L o n t round ’I' '• o :: • ■ e r i. e n t v; i t h L ii e c u o f f i c i c n t 



d e V e 1 o p o d fro m 


the simulation. 


t !i 


e lattice sp 


a c i n g 


0 f 


L J e 


base material. 


tlie temDerature 


o f 


the s y s t e m , 


t h G r 


:: a s 


s o f 


the diffusing 


particle, and 


the 


force cons 


t a n t 


of 


t h e 



harmonic oscillator approximation to the potential function 
b e t V/ e e n the diffusing atom and the base material must i; e 
provided. As an ex a'm pie, the (iiffusion coefficient for 
carbon diffusing through, fee iron at 910 degrees Celsius has 

_ •7 9 

been e x p e r i m e n t a 1. 1 y determined to be 1.91 x 10 ' c m “ / s 

(Askeland, 198 4). The lattice unit for fee iron is 3.580 x 

_ O 

10 cm (Askeland, 19 S 4) , the force constant for carbon- 
iron is in the range 2.6 5 x 10^ to 17.73 x lO’’ dynes/cm 
( A . J . Gordon, 1 972), and the atomic v; e i g h t of carbon is 
12.011 gr/nol. Frora these constants, the experimentally 
determined diffusion coefficient listed above is 
transformed into a range of values that fall b e t e e n 0.23 
and 0.61 i n t e r - p 1 a n a r distance squared per t i m e s t e p . The 
diffusion coefficient, values deter r; ined by t’ne time 
averaging over 100 t i in e s t e p s , fall between 0.19 and 0.83 
inter - planar distance squared per timestep. Thi(^ values for 
the c o e f f i c i. e n t d e t e r m i n e d b y the s i m u 1 a t i o n ’ s t i m e 
a v e r a g L II g over 1 0 C 0 t i m, e s t e p s is between 0.2 6 a n d 0.42 
i n t e r - p 1 a n a r distance squared per timestep. This indicates 
that tine simulation is determining the diffusion coefficient 
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within acceptable limits. 



3 . Concentration r o £ 1 1 C 

i 1 G concentration p r o f i ] e across a :r. e .t b r a n c s C oi' Iv- 
or o c e e d t o a linear slope across the p 1. a n e s of t h. e c r y s t 1 
structure. The s i ';■! u la t i o ri 'DIFrUSE' be 3 ins to aporoacl. t'.iO 
correct forr. at th. e 1-1000 tine average of the concentre ti. on 
(see figure 14). It then proceeds to oscillate about th. e 
infinite tine profile result. (see figure 15-17). Thi". 
oscillation about the infinite tine line indicates that ui:e 
simulation reaches an approximate steady-state 
configuration. To test this T. y p o t h e s i s , the 1 9 0 0 1 - 2 0 C) 0 0 
time average profile was statistically ch. ecked against tiie 
linear infinite fit. The statistical method useci v;as a 
linear regression of the data to the infinite result. T'*. e 
regressionre suits gave a probability of 0.999 that tiie data 
does fit tlie linear slope of the infinite time profile. Th. e 
sane approximate r e s u 1. 1 was found for the 1 4 9 0 1 - 1 5 C 0 0 and 
the 6 0901-61000 tine average concentration profiles. ThvO 
same regression was fit, using the data averaged over 
2 0 0 0 0 t i m e s t e p s , found that the probability that its 
concentration profile could bo fit to t!ie infinite ti'ic 
profile was 0.953. 

because of t:.e concentration profile's closeness of 
fit to the infinite time profile, the apparent derivation 
of the actual diffusion coefficient by the s i n; u 1 a L i o n , a n c. 
tlio simulation's ability to obey tiie conservation of matter. 
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Figure 14. 
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Figure 15. 
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Figure 17. 
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it is roasonablG to assert that the simulation is acting; as 
as an a p p r o i t: a t e physical mode] of the interstitial 
diffusion mechanism, throu.eh a fee crystal. 



CALCULATIOi: OF THE DIFrUSION COEmCIEE'T 



1 , F i c ic s First L a liethod 

'DIFFUSE’ initially develops the diffusion 
coefficient by usin^ the Pick’s first law method. Inis 
n e t hi o d t a 1: e s the time average flux and divides it by the 
time average first derivative of the concentration, iiy this 
method, th. e diffusion coefficient shows \v’hat appears to be; 
an oscillation across the planes (see Figures 13 c. 19). The 
cause of these oscillations could be th.e result of either 
noise in the simulation, the data fitting a horizontal line 
across the plane, or* a positional dependence on the 
diffusion coefficient. To test whether th. e data fit a 
horizontal line across the plane, a linear regression was 
fit to the data. The results of the fit are plo.tted in 
Figures 13 T 19. The linear correlation factor for the time 
a V e r la g e gives a p r o b a b i 1 1 i t y of 20 - 30 per cent t hi a t t ■ : e 

data fits a line of constant slope and less then 2 0 
percent probability that the data fits a horizontal line 



across the planes. 

To investigate the chance that the oscillation was 
noise and not a positional dependence on th. e diffusion 
coefficient, the diffusion coeficient data as a function of 
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DIFFUSION COEFFICIENT 
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r i. g u r e 19. 



50 



PLANES 



1 1: is apparent iron: the t o t a ’n 1 g s that t h o results of l ■. e 
L o "1 e t li o s are e r a s t i c a 1 1 y c‘ 1 f f a r e n t . h n 1 i r o t 't £' ; r i 

r.'. e t i . o d data 1 1: e second 1 a data is not equally d i s t r i I.. u ted 
a b o u ti t s ;n e a n and all the data are not v; i t h i n t c s t a n dare 
deviations of the mean. The means of the t o r.: e t h o d s are 
also not v; i t h i n two standard deviations of each other. 
Therefore, the diffusion coefficient nay be dependent on the 
position of tlie diffusing atom. TIi e next stop in Llio test 
of the hypothesis is to find t !i e first derivative of t h. c 
diffusion coefficient and compare the time derivative of ti. e 
concentration and the result of Pick’s second lav; \;ith tiio 
)i on - linear diffusion coefficient (Equation 3a). The results 
of this comparison can be seen in Tables S and 9. 



TAF/lE 8. 

TAFLE OF SHORT TliiE-A VERAGED TIME DERIVATIVES 



Plane 


T i rr. 0 Derivative 


Pick’s S e c o 


U 


.00053 


- .02396 


/ 


- .00056 


-.00320 


-j 


- . 00222 


- .00199 


0 


-.00111 


-.00481 


7 


.00000 


- .00173 
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i.O'jT. : Table of tine derivatives from 

b V ’ n u a t i o n 3 a u s i n q the a. v e r a r> o s 
t i n o s L e p s . 



s i r; u 1 a L i o n a n d (’. CT i e d 
over t h c ] 9 9 0 1 - 3 E‘ 8 ' ’ T 



3 LE 



o 

✓ 



ta:lz of lo;:g ti..e-av!:fagef ti; 



_..i A i i X 



Plane 

a 

4 

5 

6 

7 



L r': o D e r i V a L i V e 
.COO 17 
.00000 
-.00006 
-. 0001-1 
-.00011 



" i c ' s S e c o n L L ; 



-.02630 

- .00142 

- .000 2 7 
.00104 

-.00107 



!.'0TF: Table of time derivatives from simulation and derived 

hy Er nation 3 a using the averages over tiie ] 0 0 0 1 - 2 O' 0 0 
tiniestens. 



T b. o results of this comparison s h o v; that 1 1: o non-linear form 
of Fic!:*s second lav; does not agree with the simulation 
results of tlie tine derivative of the concentration. This 
•suggests that Pick’s first lav;, taken literally, nay not 
be the best method to fit the s i m. u 1 a t .1 (J n data to a diffusion 
coefficient. It may be the diffusion coefficient sriould be 
v; i t h i n the partial derivative. 

3 . F o k k e r - ? 1 a n !: k e t h o d 

The last method used to s e a r c li for a positional 
dependence v; a s to attempt to fit a function t. o the d i f f u s i o n 
coefficient and solve the F o k k e r - P 1 a n equation ( F c u a t i o n 
13) v; it: this function and t 'n o results of the s i : ^ i: 1 a t i c n . 

T ii e folio 1. n g functional for n v; a s assumed, b nsec! on t h 
d i f f u s i o n coefficient resulting from t h. e 1' i. c i; ’ s first la w 
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met hod : 



!)(::, t) = A + i;sin(Cx + Jt + A) 



( 12 ) 



Using tr..is ecuation and Llie resulLs fror. tho s i r. u 1. a I: i o , L'lc 
To ':ker-Plan;c equation ta!ces the follo’.;ing for:.’.: 



9C/9t = Bsin(Cx + Dt + E)[0.5C^(CVOL) + C(:'ERCl) = .5(Li:r.C2)] + A/2 

(13) 



The I'ASL library subroutine ZXS50 \:l\s used to evaluate th.e 
Folcker -Plank equation. ZXSSC finds the ninini un of tlic sun 
of squares of M functions in X variables using a finite 
difference Levenberg-Marquardt algorithm. In the simulation 
the N variables are tlie A, B, C, D, and E variables of th. c 
diffusion coefficient function, and the A functions are 
Equation 13 evaluated at different timesteps fror: tiic 

simulation. 

It v/as determined that the best fit for the 
coefficients A , B , C , D and E was w 'i e n t b. o y all equaled 
zero. This result indicated that the oscillations in the 
diffusion coefficient ’./ere probably based on the numerical 



noise in the s v s t e n . 



n o 



clear 



V. co::cLUSio;;s a'.’d ?.i-:co:;:a:::'j;.TTo;;: 



T' h e (i a t r. from the simulation 5 a v e 
determination, one v/ay or tlie other, about the positional 
dependence of the diffusion coefficient. The sim. ulat.ion 
should be further improved to malce a specific deter min a ti. on 
of the positional de'n end once of the diffusion coefficient. 
The follov;inp, improvements to the simulation sliould give 
it the ability to make the determination of the dependence. 

First, the time averaging of tl-ie flux, concentration, 
ant; derivatives nay have averaged out any real positional, 
dependence. Therefore, the niimber of timesteps thiO 
siraulation averages over should be shortened, or deleted 
altogether. Instead, the few times that the derivatives are 
zero should be ignored and the diffusion coefficient should 
be inspected at all other times and planes. 

Second, the size of the array should be increased. Th.o 
size used in this thesis, 6 X 6 X 16 , may not have been 
large enough. It appears that v; i t h this size a single 
p a r t i. c 1 e change in a plane loads to large changes in the 
nar a meters determined by the simulation. The array should bo 
increased in size to at least a 14 X 14 X 5 0 a r r a y . This 
size should na!:o a single particle ciiango in a plane equal a 
o n e percent change in the determined factors of 1 1 > o 



J \j 



simulation. 



Third, t ’n e functional f o r r.: of L !i c ri i f f u s i o n. 

coefficient, and tie initial values for L h: e analysis of L ’ e 
unkov;n variables for the F o Ic Ic e r - ? 1 a n :: equation need to be 
investinated furtlier. 

A major result of this thesis is its ir. provement in the 
bookeepin.q ability of diffusion simulations. 'DIFFjSF' uses 
aTrays to keep track of values that are needed to simulate a 
i'.i o r e s o p h i s t a c a t e cl model of diffusion. These arrays k, e c p 
trade of t ii e n e a r e s t - n e i 3 h b o r s for the diffusing atoms, and 
t 'n e planes that the atoms are in. T ri i s a 1 1 o ’>/ s t o 
simulation to be C'/ritten more simply than earlier 
s i !u u 1 a t i o n s and a 1 1 o v; s t e more sophisticated n o d o 1 i n of 
the face-centered-cubic crystal. s to be simulated. 

The simulation's improvements to the method by \;hicli tb. o 
diffusion process is modeled, along v; i t h these s u ", g c s t e d 
improvements to the simulation, should a 1 1 o v; t ii e simulation 
to determine if the diffusion coefficient has a positional 



d o T' e n d e n c e 



The f o 1 1 o \v i n 2 



names is set u n to 



appt;;dp:: a 

V APT ABLE GLOSS A BY 
glossary of variable 
aid in using the computer programs ' D I F S E T ’ and ' D I F F j S . 
The glossary consists of three columns, the nar.e of the 
variable, the type of variable and t li e physical m c a n i. n g of 
the variable. 
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VARiAr;Li: 


i . n . i 




VAX I ALL';: type 


PEYSICAL year: 


Acc::s 






2-D A PE AY 


AVER AG E 










co:;cE::Tr:AYioR 
EACE PEAEE. 


ACVOL 






2-D ARRAY 


AVERAGE 
COECEETRATIOE 
PER UEIT VOL. 

r A r !] pT/'-- 

Ij ^ 1 1 1 i j / t ^ . Lj • 


ADEPXl 




• 


2-D ARRAY 


A. V E R A G E 1 

D E R I V A. T I V E. 
TEE CCEiCEEi: 
T T 0 A T E A 

RLAEE. 


ADERC2 






2-D ARRAY 


A V E R A G E 2 

D E R I V A T I V E 
EACH PEAEE. 


ADEPT 






2-D AP.PAY 


A V r. R A G E T I 
D r. R I V ATI V E 

i . . i j w V . * . * i i 

T T 0 E A T E A 
PEAEE. 


A FLUX 






2-D ARRAY 


AVERAGE ELEE 
EA'CE PEAEE. 


AT:iSTP 






ARGUIIEMT 


USED A S 1 
ARGUETET IE 
THE A V E R A G 
ARRAYS . 


ATI 






PARAMETER 


E U E R E R. 
TIEESTEPS r 
T El F, A V E P. G 
ARRAYS . 


AVGTU 






PARAi.'ETER 


i: U :i B E R 
TIEiESTEPS iC 
AVERAGED OVER 


avgt::i 






PARAMETER 


PARAEETEP IE' 
TO EiOLD EAST 
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VA"! IATj-L" 



va?.ia:l:: ty? 



P'lYSiCAL 



C^.A' 



0 : ■ > •' r) \ 

^ j 1 a . i . ; i i 



co:;cF.:.Tr.ATio.: 

FOF FAc:; pla: F 

AFj TIF'ESTFF. 



CVOL 



^ - J.' n P. .-I 1 



cc::ce::tfat:o:; 

PFF U::iT VOL AT 
FACT. PLATE ATP 
TITESTFP. 



1 V . . * . i J » . 



O n \ -n -n \ y 

Z. — IJ A { w iv 1 



T’EAPPFT-TFIGTrOF 
TA 2 LF . 



CFYST'^ 



O T-, \^T \ 

^ — u A 1 V . . V J 



TTI'ET JITT :FI 0 T- 
AL ARE AY OF FCC 
CRYSTAL . 



CRY ST 



1 -D ARRAY 



o;;r-di:;ftsio:;al 

EOillVALFTT TC 
CRYST 3 . 



r n r 1 



^ - iy A .V n i 



1 ST jERIVATIVF 
OF TEE CCTCFT- 
T RAT TOT AT FACT 
PLATE. 



DFRC 2 



o F', r> T? A 

^ L) 1 k i V iV a 1 



2 TD BFRIVATIV 
AT EACi: PLATE. 



n r D T 



2 -D ARRAY 



T I E D E I V i\ T I E 
AT EACH PLATE. 



JITCOI 



2 -D ARRAY 



DIFFUSTCT SOP;.-- 
FICTLTT at E.IC.- 
? L A T T , L E T E P - 
TTTE- RY FICT’T 
FIRST LAI . 



TY 



IFC 02 



O 'T P \ Y 



DIFFUSiOT COFF- 
FICIi'FT RETER- 
TTTED RY F ICE'S 
SECOND LA'.;. 



FLTT 



2 -]-: A ’ lAY 



nux acr:;ss 

PLATE AER T 
TE.P. 



FACT 
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CT) 







eaxeye:^ 


■•p-.:er (r: i a : 

T 1 0 P . 


1 . . 1 




PAEA.X'Ti-.E 


.. b . -'L.. u.‘ _ E_:v 
TE.AR TPAE. 


IS "EDI 




D ' “D \ ’ : r *'r t: p 


R a::i)0p ::u ' r :? 

C ” T? p 

kJ 1 j i j U • 


ISEED2 




V p ' t f n 'Y' r 


RAEDOP EPPEER 
SERD. 


J'.IAX 




PAEAXETEE 


PUPRER OP-PLAERS 
IE' Y DIREC'i’KPP 


J'.iAXl 




PAEAMETEE 


OEE LESS TEAE 
JPAE. 


T'’ 

1 «. 




argu;:e::t 


ARGLPIEET TC' 'LL 
2-DIP'EESIOEAL 
A R R A Y s . n E. P : " - 
SEEiTS TE>E PLAE", . 


r ' : A V 

L\ L i i 1 < L 




PAEAPETER 


EUP5EE or PLA'fEE 
IK TEE Z riREC- 

1 L . . • 






PARAPETER 


Ki ^'v • 


L 




ARCPJPTPT 


ARGPPE"T E E 
ARRAYS ZPLE A P 
CIPIST . 


l;:ax 




PA RAPE TER 


::upp.!:r op stv"s 

IE CRYST. 


IXiAX 1 
E..".: 




PARAPETER 


LPAE-ET EPAE . 

GEE PLA"E p;:rs 
TEAE’ ALL T' E 
SITES IE (P’VS';. 




ARGLP.EE'i 


\ ^ r ] T ' ■ ■' { ’ ' ’ ■ ^ 

4 ^ ‘ j L . . j . . 1 



01 



L;l 



V;. niAi'.!.!". 



.)T 



iVSICAl. 



LI 



lavgt:! 
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CUTrUT 



n ;■ 

iV 1 






no:: 2 



• SiP 



T 1 



;?l:i 



.1 . Li L . 



P \V \ 
1 iv n 



ETEH 



p A \ * p T 



LOGICAL 



1-D APvRAY 



1 r. A 13 A 
i “ t\ L\ i\ V I 1 



.r.GLL'iE'.IT 



PA 



;ahet 



■ R 



l-D A 



RAY 



0 2 



1 . . r w i-. .. 1 j . 

' • ~ - T* T“ ' 

. » — I i . 1 _ L. i . . » ) 

TO ri;:o 
:Tic:r;''Oi: . 



LAr.(;LL i o .. . 

OF Tin: average: 
TIEESTEFS . 



LARGEST EELIRi" 
OF TIEESTEi'S. 

o'ETErEIEES 

a:iole:t of outprt 

TEE ? R 0 ' - E A . . 
GIVES . 



ARRAY OF EAERO: 
lEJEEERS USED FO: 
PICE I EG CEYs: 
SITE TO R : 
LOOSED AT. 



ARRAY OF R A EDGE 
EUEDERS USED TO 

p I c i: :: e a r r. s t - 

LEIGH DOR TO EE 

look::d at. 



AR G U E E E 



r ■ 

i V./ 



r T • ' V T. r' ^ ^ 1 f • - •' r O 

i- ^ j ^ I j i. j ■ j> ^ 

AE'D FEET. 

FIRS T T 1 EES 

V r t: r 

progra;; . 



f* D 



U 1 



»li v/t. 



VS 



r>T \ ■ f C 

*. Lj c \ . 1 




LISTi;;G 



A?P":;jix 
'DIFSET' CODE 
TIic folio v;ing is t'ne computer code i 
program ' D I F S E T ' . This code develops the 

program 'DIFFUSE'. 



n 

J 



or t a c s 0 t - 
input for t 



C 'DIFSET* 

C MARK R, POLNASZEK 

C JUNE 1986 

C 
C 

C •DIFSET* HAS DESIGNED TO INTIALIZE THREE ARRAYS THAT ARE NEEDED 
C IN THE PROGRAM 'DIFFUSE* • THESE ARRAYS CRYST(L), CRNNBR(L,R), 

C AND ZPLN(L) ARE USED IN 'DIFFUSE* TO SPEED UP THE CALCULATIONS 
C OF THE DIFFUSION PROCESS. 

C 

C 

C THE FOLLOWING DEFINES THE VARIABLES OF THE PROGRAM AS REAL OR 
C INTEGER VALUES, AND DIMENSIONS THE ARRAYS* MEMORY SIZE. 

C 

C 

00020 INTEGER PLNMAX,CRYST( 2000 ) , CRYST3I 10 , 10 ,20 ) ,CRNNBR( 2000 ,12 ) ,R, 

1 ZPLN(2000), ISEEDl, ISEED2, FLUX(2000), CONS ( 2000 ) ,MTMSTP , 

2 T1 , AVGTM , PLSTM , OUTPUT 



C 

C 

C THE FOLLOWING READS IN THE INPUT TO THE PROGRAM AND ECHOS IT 
C BACK TO AN OUTPUT FILE. 

C 

C 

READ (5, *)IMAX, UMAX, KMAX,MTMSTP ,T1 , ISEEDl ,ISEED2 , AVGTM , PLSTM, 
1 OUTPUT 

WRITE J6,26) IMAX, UMAX, KMAX,MTMSTP ,T1 , AVGTM 
00026 FORMAT ( *1* ,618) 

WRITE 16,28) ISEEDl, ISEED2 
00028 FORMAT 1*1* ,2116) 

C 

C 

C THE CODE BELOW INITIALIZES CONSTANTS THAT ARE USED THROUGH OUT 
C THE PROGRAM. 

C 

C 

00040 PLNMAX = IMAX * UMAX 
00050 LMAX = KMAX * PLNMAX 

PRINT 52, PLNMAX, LMAX 
00052 FORMAT I *1* ,218) 

C 

C 

C THE FOLOWING DO LOOPS SET UP THE FIRST ARRAY ZPLN(L). THIS ARRAY 

C HOLDS THE VALUE OF THE PLANE NUMBER FOR THE GIVEN SITE NUMBER, L. 

C 

C 

DO 3192 L = 1,LMAX 
ZPLN(L) = (INT (L/PLNMAX)) + 1 
DO 3191 N = PLNMAX, LMAX, PLNMAX 
ZPLN(N)= N/PLNMAX 

03191 CONTINUE 

03192 CONTINUE 
C 

C 

C THE FOLOWING DO LOOPS ARE USED TO SET UP THE THREE DIMENSIONAL 
C ARRAY CRYST3(I,J,K). THIS ARRAY HOLDS THE VALUE OF THE TYPE 
C OF ATOM IN THE CRYSTAL SITE OR THE FACT THE SITE IS EMPTY. 

C 

C 

00060 DO 102 K = 2, KMAX, 2 
00062 DO 70 J = 1, UMAX, 2 
00064 DO 68 I = 1, IMAX, 2 

00066 CRYST3(I,J,K) = 0 

00068 COr4TINUE 

00070 CONTINUE 

00072 DO 80 J = 2, UMAX, 2 
00074 DO 78 I = 2, IMAX, 2 



00076 

00078 

00080 

00082 

00084 

00066 

00088 

00090 

00092 

00094 

00096 

00098 

00100 

00102 

00104 

00106 

00108 

00110 

00112 

00114 

00116 

00118 

00120 

00122 

00123 

00124 

00125 

00126 
00128 
00130 
00132 
00134 
00136 
00138 
00140 
00142 
00144 
00146 
00148 
00150 
00152 
00154 
00156 
00158 
00160 
00162 
00164 
00170 
00172 
C 

C 
C 
C 
C 
C 
C 
C 

00174 

00176 

00178 

00180 

00182 

00184 

00186 

00188 

00210 

C 

C 



CRYST3(I,J,K) = 0 
CONTINUE 
CONTINUE 

DO 90 J = 2, UMAX, 2 
DO. 88 I = 1,IMAX,2 
CRYST3(I,J,K) = 2 
CONTINUE 
CONTINUE 

DO 100 J = 1, UMAX, 2 
DO 98 I = 2,IMAX,2 
CRYST3(I,J,K) = 2 
CO^mNUE 
CONTINUE 
CONTINUE 

DO 124 K = 1,KMAX,2 
DO 114 J = 1,JMAX,2 
DO 112 I = 1,IMAX,2 
CRYST3(I,J,K) = 2 
CONTINUE 
CONTINUE 

DO 123 J = 2, UMAX, 2 
DO 122 I = 2,IMAX,2 
CRYST3(I,J,KJ = 2 
CONTINUE 
CONTINUE 
CONTINUE 

DO 146 K = 3,KMAX,2 
DO 134 J = 2, UMAX, 2 
DO 132 I = 1,IMAX,2 
CRYST3(I,J,K) = 0 
CONTINUE 
CONTINUE 

DO 144 J = 1, UMAX, 2 
DO 142 I = 2,IMAX,2 
CRYST3(I>J,K) = 0 
CONTINUE 
CONTINUE 
CONTINUE 
K = 1 

DO 158 J = 2, UMAX, 2 
DO 156 I = 1,IMAX,2 
CRYST3(I,J,K) = 1 
CONTINUE 
CONTINUE 

DO 172 J = 1, UMAX, 2 
DO 170 I = 2,IMAX,2 
CRYST3(I,J,K) = 1 
CONTINUE 
CONTINUE 



END OF THREE DIMENSIONAL CRYSTAL SET UP AND THE BEGINNING OF 
DIMENSIONAL EQUIVELENCE. THE ONE DIMENSIONAL ARRAY CRYST( L ) 
ARE IDENTICAL TO THE THREE DIMENSIONAL CRYST3( I , J,K ) . 



DO 188 K = 1,KMAX 
DO 186 J = 1,JMAX 
DO 184 I = 1,IMAX 

L = I ♦nj-11^IMAX) + ( (K-1)*PLNMAX) 
CRYST(L) =CRYST3(I,J,K) 

CONTINUE 
CONTINUE 
CONTINUE 
KMAXl = KMAX -1 



THE ONE 
VALUES 
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C THE FOLLOHING GROUP OF DO LOOPS DEVELOP THE TWO DIMENSIONAL 
C ARRAY CRNNBR(L,R). THIS ARRAYS VALUES ARE THE TWELVE NEAREST 
C NEIGHBORS TO THE CRYSTAL SITE L. 

C 

C 



00220 


DO 380 K = 2,KMAX1 




00230 


J = 1 




00240 


11 

M 




00250 


L = I ♦ 


(K-1)»PLNMAX 


00260 


CRNNBR(L,1) = 


IMAX ♦ K*PLNMAX 


00270 


CRNNBR(L,2) = 


I ♦ (JMAX-maMAX ♦ K^PLNMAX 


00280 


CRNNBR(L,3) = 


I+l ♦ K»PLNMAX 


00290 


CRNNBR(L,4) = 


I ♦ IMAX ♦ KJfPLNMAX 


00300 


CRNNBR(L,5) = 


IMAX ♦ ( JMAX-I)*IMAX ♦ (K-I)»PLNMAX 


00310 


CRNNBR(L,6) = 


I + l ♦ (JMAX-H»IMAX ♦ (K-I)»PLNMAX 


00320 


CRNNBR(L,7) = 


I + l ♦ IMAX ♦(K-I )*PLNMAX 


00330 


CRNNBR(L>8) = 


IMAX ♦ IMAX ♦ (K-I)»PLNMAX 


00340 


CRNNBR(L,9) = 


IMAX ♦ (K-2)»PLNMAX 


00350 


CRNNBR(LaO) 


= I ♦ ( JMAX-I )*IMAX ♦ (K-2)»PLNMAX 


00360 


CRNNBR( L,ll) 


= 1*1 ♦ (K-2>»PLNMAX 


00370 


CRNNBR(l/12) 


= I ♦ IMAX + (K-2)*PLNMAX 


00380 


CONTINUE 




00390 


DO 530 K = 2,KMAX1 




00400 


J = 1 




00410 


I = IMAX 




00420 


L = I ♦ (K-1)*PLNMAX 


00430 


CRr^r^BRtL,!) 


= I-l ♦ K^PLNMAX 


00440 


CRNNBR( LfZ) 


= I ♦ ( JMAX-1)*IMAX ♦ K«PLNMAX 


00450 


CRNNBR( L,3) 


= 1 ♦ K*PLNMAX 


00460 


CRNNBR(L,4) 


5 T ♦ IMAX +K*PLNMAX 


00470 


CRNNBRl L,5) 


= I-l ♦ ( JMAX-1 )*IMAX ♦ (K-1 )^(PLNMAX 


00480 


CRNNBR( L,6) 


= 1 ♦ ( JMAX-1 )*IMAX ♦ (K-1)*PLNMAX 


00482 


CRNNBR( L,7) 


= 1 ♦ IMAX ♦ (K-1)*PLNMAX 


00484 


CRNNBR( L,8) 


= I-l ♦ IMAX ♦ (K-1)*PLNMAX 


00490 


CRNNBR( L,9) 


= I-l ♦ (K-2)*PLNMAX 


00500 


CRNNBR( L,10) 


= I ♦ ( JMAX-1 )^IMAX ♦ (K-2)*PLNMAX 


00510 


CRNNBR( L,ll ) 


= 1 ♦ (K-2)*PLNMAX 


00520 


CRNNBR( L,12) 


= I ♦ IMAX ♦ (K-2)*PLNMAX 


00530 


CONTINUE 




00540 


DO 700 K = 2,KMAX1 




00550 


J = UMAX 




00560 


I = IMAX 




00570 


L = K* PLNMAX 




00580 


CRNNBR(L,1) 


= I-l ♦ (J-1J*IMAX ♦ K^PLNMAX 


00590 


CRNNBR( L,2 ) 


= I ♦ (J-2)*IMAX ♦ K*PLNMAX 


00600 


CRhJMBRl L,3 ) 


= 1 ♦ (J-1)*IMAX ♦ K^PLNMAX 


00610 


CRNNBR( L,4) 


= I ♦ K*PLNMAX 


00620 


CRNNBR( L,5) 


= I-l ♦ (J-2»*IMAX ♦ (K-1)*PLNMAX 


00630 


CRNNBRl L,6 ) 


= 1 ♦ (J-2)*IMAX ♦ (K-1)*PLNMAX 


00640 


CRNrJBR( L,7) 


= 1 ♦ (K-1J*PLNMAX 


00650 


CRNrBRt L,8) 


= I-l ♦ (K-1J*PLNMAX 


00660 


CRNNBR(L,9) 


= I-l ♦ (J-1J*IMAX ♦ (K-2)*PLNMAX 


00670 


CRNNBRt L»10 ) 


= I ♦ (J-2)*IMAX ♦ (K-2)*PLNMAX 


00680 


CRNNBRt L,ll ) 


= 1 ♦ (J-1)*IMAX ♦ (K-2J*PLNMAX 


00690 


CRNNBR( L,12) 


= I ♦(K-2)*PLNMAX 


00700 


CONTINUE 




00710 


DO 850 K = 2,KMAX1 




00720 


J = JMAX 




00730 


1 = 1 




00740 


L = I ♦ (J-1)*IMAX ♦ (K-1)*PLNMAX 


00750 


CRNTBRl L,1 ) 


= IMAX ♦ (J-D^^IMAX ♦ K^PLNMAX 


00760 


CRNriBR( L,2 ) 


= I ♦ (J-2»*IMAX ♦ K^PLNMAX 


00770 


CRNNBRl L,3 ) 


= I-H ♦ (J-1)*IMAX ♦ K^PLNMAX 


00780 


CRNNBR( L,4) 


= I ♦ K*PLNMAX 


00790 


CRr;NBR( L,5 ) 


= IMAX ♦ (J-2J*IMAX ♦ (K-1 )^^PLNMAX 


00800 


CRNtJBR(L,6) 


= I+l +(J-2)^IMAX ♦ (K-1)*PLNMAX 


00810 


CRNNBR(L,7) 


= I + l +(K-1 )^^PLNMAX 


00820 


CRNNBR( L,8) 


= IMAX ♦ (K-1)*PLNMAX 
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00830 

00840 

00842 

00844 

008S0 

00860 

00870 

00880 

00890 

00900 

00910 

00920 

00930 

00940 

00950 

00960 

00970 

00980 

00990 

01000 

01010 

01020 

01030 

01040 

01050 

01060 

01070 

01080 

01090 

01100 

OHIO 

01120 

01130 

01140 

0x150 

01160 

01170 

01180 

01190 

01200 

01210 

01220 

01230 

01240 

01250 

01260 

01270 

01280 

01290 

01300 

01310 

01320 

01330 

01340 

01350 

01360 

01370 

01380 

01390 

01400 

01410 

01420 

01430 

01440 

01450 

01460 

01470 

01480 



(K-2 I^PLNMAX 
(K-2 I^PLNMAX 
(K-2)^^PLNMAX 



CRNf^R(L,9) = IMAX ♦ (J-1)*IMAX ♦ 
CRNNBR(L,10) = I + (J-2)* IMAX ♦ 
CRNNBR(L,11) = I+l ♦ (J-l)^IMAX ♦ 
CRNNBR(L>12) = I + (K-2)*PLNMAX 



CONTI r>»UE 
JMAXl = JMAX -1 
IMAXl = IMAX -1 

DO 1050 K = 2,KMAX1 
1 = 1 

DO 1040 J = 2, JMAXl 

L = I ♦ (J-l)^IMAX + (K-1)*PLNMAX 

CRNNBR(L,1) = IMAX + (J-1)*IMAX + tK^^PLNMAX) 
CRNNBR(L,2) = I + (J-2)*IMAX + K*PLNMAX 
CRNNBR(L,3) = I + l + ( J-1 )^IMAX ♦ K^PLNMAX 
CRNNBR(L,4) = I + J^^IMAX ♦ K^PLNMAX 
CRNNBR(L,5) = IMAX + ( J-2 )*IMAX ♦ (K-l)^PLNMAX 
CRNt4BR(L,6) = I + l + (J-2)^^IMAX ♦ (K-D^PLNMAX 
CRNraR(L,7) = I + l +(J*IMAX) + (K-1)*PLNMAX 
CRNNBR(L,8) = IMAX + J^^IMAX + (K-D^PLNMAX 
CRNNBR(L,9) = IMAX + (J-D*IMAX + (K-2)^PLNMAX 
CRN^4BR(I,10) = I + (J-2)*IMAX + (K-2)*PLNMAX 
CRNNBR(L,11) = I+l + (J-l)^IMAX + (K-2)*PLNMAX 
CRNNBR(L,12) = I + J^IMAX + (K-2)^PLNMAX 



CONTINUE 

CONTINUE 

DO 1230 K = 2,KMAX1 
J = 1 

DO 1220 I = 2, IMAXl 

L = I + (K-1)*PLNMAX 
CRNNBR( L,1 ) = I-l 
CRNNBRt L,2 ) = I + 
CRNNBRt L,3 ) = I + l 
CRNNBR(L,4) = I + 
CRNNBR(L,5) = I-l 
CRNNBRtL,6) =1+1 
CRNt4BR(L,7) = I + l 
CRNNBRt L,8) = I-l 
CRNNBR( L,9) = I-l 
CRNNBR(L,10) = I 
CRNNBRIL,!!) = 1+ 
CRNNBR(L,12) = I 

CONTINUE 

CONTINUE 

00 1410 K = 2,KMAX1 
I = IMAX 

00 1400 J = 2, JMAXl 

L = I + IJ-1)*IMAX + 
CRNNBR( L,1 ) = I-l 
CRNNBR( L,2 ) = I + 
CRNNBR(L,3) = 1 + 
CRNNBR( L,4 ) = I + 
CRNNBR( L,5) = I-l 
CRNNBR( L,6 ) = 1 + 
CRNNBR( L,7) = 1 + 
CRNNBR( L,8) = I-l 
CRNNBRt L,9) = I-l 
CRNNBR(L,10) = I 
CRNNBRI L,ll ) = 1 
CRNNBR( L,12 ) = I 
CONTINUE 
CONTINUE 

DO 1590 K = 2, KMAXl 
J = JMAX 

DO 1580 I = 2, IMAXl 

L = I + ( J-1 I^IMAX + 
CRNNBR( L,1 ) = I-l 
CRNNBRt L, 2 ) = I + 
CRNNBRt L,3 ) = I+l 



+ K*PLNMAX 

tJMAX-imiMAX + K*PLNMAX 
+ K*PLNMAX 
IMAX +K^^PLNMAX 

+ tJMAX-D^IMAX + tK-l)^^PLNMAX 
+ t JMAX-1)*IMAX + tK-l)^^PLNMAX 
+ IMAX + tK-ll*PLNMAX 
+ IMAX + tK-l)*PLNMAX 
+ tK-2 )^^PLNMAX 

tJMAX-D^IMAX + tK-2 I^^PLNMAX 
+ tK-2)*PLNMAX 
IMAX + tK-2)*PLNMAX 



tK-1 )^PLNMAX 
+ tJ-l)^IMAX + K^PLNMAX 
tJ-2)*IMAX ♦ tK*PLNMAX) 
tJ-l)*IMAX ♦ K*PLNMAX 
J^^IMAX) ♦ K^PLNMAX 
+ tJ-2)*IMAX + tK-l)*PLNMAX 
tJ-2)*IMAX + tK-ll^PLNMAX 
J^^IMAX + tK-l»*PLNMAX 
+tJ*IMAXI ♦ tK-ll*PLNMAX 
+ tJ-l)+IMAX + tK-2)*PLNMAX 
tJ-2)*IMAX ♦ tK-2 )*PLffMAX 
tJ-l)*IMAX ♦ tK-2)^PLNMAX 
J^IMAX ♦ tK-2l^PLNMAX 



tK-1 )*PLNMAX 
+ tJ-ll^IMAX +K*PLNMAX 
tJ-2)*IMAX + K^PLmAX 
+ tJ-l)*IMAX + K^PLNMAX 



01^^90 

01500 

01510 

01520 

01530 

01540 

01550 

01560 

01570 

01580 

01590 

01600 

01610 

01620 

01630 

01640 

01650 

01660 

01670 

01680 

01690 

01700 

01710 

01720 

01730 

.01740 

01750 

01760 

01770 

01780 

01790 

01800 

01810 

01820 

01830 

01840 

01850 

01860 

01870 

01880 

01890 

01900 

01910 

01920 

01930 

01940 

01950 

01960 

01970 

01980 

01990 

02000 

02010 

02020 

02030 

02040 

02050 

02060 

02070 

02080 

02090 

02100 

02110 

02120 

02130 

02140 

02150 

02160 



CRNNBR(L,4) = I ♦ 
CRNNBR(L,5) = I-l 
CRNNBR( L>6 ) = !♦! 
CRNNBR(L,7) =1+1 
CRNNBR(L,8) = I-l 
CRNNBR(L,9) =1-1 
CRNNBR( L,10 ) = I 
CRNNBRC L,ll ) = 1 + 
CRNNBRC L,12 ) = I 
COffTINUE 
CONTINUE 

DO 1780 K = 2,KMAX1 
DO 1770 J = 2,JMAX1 
DO 1760 I = 2,IMAX1 
L = I + (J-1)*IMAX ♦ 
CRNNBR(L,1) = I-l 
CRNNBRC L,2) = I ♦ 
CRNNBRC L,3 ) = I + l 
CRNNBR(L,4) = I ♦ 
CRNNBR(L,5) = I-l 
CRNr4BR(L,6) =1+1 
CRUmRiLf?) =1+1 
CRr>ir>lBR(L,8) = I-l 
CRNNBR(L,9) = I-l 
CRNNBRC L,10 ) = I 
CRNNBR(L,11) = 1+ 
CRNNBR(L,12) = I 
CONTINUE 
CONTINUE 
CONTINUE 

DO 1960 J = 2,JMAX1 
K = 1 

DO 1950 I = 2,IMAX1 
L = I + (J-lJ^^IMAX 



(K-1 )^^PLNMAX 
(K-1 )^^PLNMAX 



K*PLNMAX 
♦ (J-2)^IMAX + 

+ (J-2>^^IMAX + 

+ IK-1)^PLNMAX 
+ IK-1)*PLNMAX 
+ (J-l)i^IMAX + (K-2)*PLNMAX 
+ (J-2)*IMAX + (K-2 )^^PLNMAX 
1 + (J-1»^IMAX + (K-2 )^^PLNt1AX 
+ (K-2)*PLNMAX 



(K-1 )*PLNMAX 

♦ (J-1)*IMAX ♦ K^PLNMAX 
(J-2)*IMAX ♦ K^^PLNMAX 

+ (J-1)*IMAX + K^^PLNMAX 
J^IMAX + K^PLNMAX 
+ (J-2)*IMAX ♦ IK-D^^PLNMAX 
+ (J-2)^IMAX + (K-1 l^^PLNMAX 
+ J^IMAX ♦ (K-1)*PLNMAX 

♦ J*IMAX ♦ (K-D^PLNMAX 

♦ (J-l)^IMAX + (K-2)^^PLNMAX 
+ (J-2)*IMAX ♦ (K-2)^^PLNMAX 

1 + (J-1)*IMAX + (K-2H^PLNMAX 
♦ J^^IMAX + (K-2)^^PLNMAX 



CRNNBR( L,1 ) = 
CRNf(BR(L,2) = 
CRmBR(L,3) = 
CRNNBR(L,4) = 
CRNNBR(L,5) = 
CRNr>IBR( L,6 ) = 
CRNNBR(L,7) = 
CRr4rJBR(L,8) = 
CRNNBR(L,9) = 
CRNNBR( L,10 ) 
CRNNBR(L,11) 
CRNI>IBR( L,12 ) 
CONTINUE 
CONTINUE 

DO 2130 I = 2,IMAX1 
K = 1 
J = 1 
L = I 

CRNf4BR(L,l) = 
CRNNBR(L,2) = 
CRNNBR(L,3) = 
CRNNBR( L,4 ) = 
CRNNBR(L,5) = 
CRNNBR(L,6) = 
CRNNBR(L,7) = 
CRNNBR(L,8) = 
CRrJNBR(L,9) = 
CRrJNBR( L,10 ) 
CRNrJBR( L,ll ) 
CRNNBR(L,12) 



♦ (J-1)*IMAX ♦ PLNHAX 



I-l 
L 
L 

I ♦ (J-2)*IMAX + PLNMAX 

L 

L 

I+l + ( J-1 )*IMAX + PLNMAX 

L 

L 

= I ♦ J^^IMAX + PLNMAX 
= L 
= L 



I-l ♦ PLNMAX 

L 

L 

I ♦ ( JMAX-1 )^^IMAX ♦ PLNMAX 

L 

L 

I+l +PLNMAX 

L 

L 

= I + IMAX + PLNMAX 
= L 
= L 



CONTINUE 
DO 2300 J 
K = 1 
I = 



= 2,JMAX1 






IMAX 



02170 

02180 

02190 

02200 

02210 

02220 

02230 

022^0 

02250 

02260 

02270 

02280 

02290 

02300 

02310 

02320 

02330 

023^0 

02350 

02360 

02370 

02380 

02390 

02400 

02410 

02420 

02430 

02440 

02450 

02460 

02470 

02480 

02490 

02500 

02510 

02520 

02530 

02540 

02550 

02560 

02570 

02580 

02590 

02600 

02610 

02620 

02630 

02640 

02650 

02652 

02654 

02660 

02670 

02680 

02690 

02700 

02710 

02720 

02730 

02740 

02750 

02760 

02770 

02780 

02790 

02800 

02810 

02820 



= I ♦ (J-1)*IMAX 
CRNNBR( L,1 ) = I-l 
CRNTBRl L,2 ) = L 
CRNNBR(L,3) = L 
CRNNBR(L,4) = I 
CRNNBR(L,5) = L 
CRNNTOR(L,6) = L 
CRNNBR(L,7) = 
CRNNBR( L,8) = L 
CRNNBRl L,9) = L 
CRNNBR(L,10) = I 
CRNNBR(L,11) = L 
CRNNBR(L,12) = L 

= 2,IMAX1 

JMAX 

= I + (J-1)*IMAX 
CRmBRl L,l) 
CRNNBRl ) 
CRNTJBRl L,3 ) 
CRNNBRl L,4) 
CRNT^R( L,5) 
CRNNBRl L,6 ) 
CRNNBRl L,7 ) 
CRmBR( L,8) 
CRNNBR( L,9) 
CRNNBR(L,10) = 
CRmBRl L, 11) = L 
CRNNBRl L, 12) = L 

CONTINUE 

DO 2640 J = 2,JMAX1 
K = 1 
1 = 1 



♦ 1J-1)*IMAX ♦PLNMAX 



♦ 1J-2)*IMAX ♦ PLNMAX 



1 + 1J-1)*IMAX ♦ PLNMAX 



♦ J^^IMAX ♦ PLNMAX 



CONTINUE 
DO 2470 I 
K = 1 
J = 
L 



♦ 1J-1)*IMAX ♦ PLNMAX 



I-l 
L 
L 

I ♦ 1J-2)^IMAX ♦ PLNMAX 

L 

L 

I+l + 1J-1)^IMAX ♦ PLNMAX 

L 

L 

I ♦ PLNMAX 



L 



CONTINUE 
1=1 
J = 1 
K = 1 



I = IMAX 
L = IMAX 



= I + 1J-1)^IMAX 

CRNNBRl L,1 )=IMAX + 1J-1)*IMAX ♦ PLNMAX 
CRNNBRl L, 2) = L 
CRNNBRl L, 3) = L 

CRNh>IBRl L,4) = I ♦ 1J-2)*IMAX ♦ PLNMAX 
CRNr>IBRt L,5) = L 
CRNNBRl L,6 ) = L 

CRNNBRl L, 7) = I+l + lJ-1 )*IMAX ♦ PLNMAX 
CRNNBRl L, 8) = L 
CRNf^BRlL,9) = L 

CRNNBRl L, 10) = I + J^IMAX ♦ PLNMAX 
CRNfJBRlL,ll) = L 
CRNNBRl L, 12 ) = L 



CRNNBRl 1,1 ) 
CRmBRl 1,2 ) 
CRmBRl 1,3 ) 
CRNNBRl 1,4 ) 
CRNNBRl 1,5) 
CRNT4BR1 1,6 ) 
CRNmRi 1,7 ) 
CRNNBRl 1,8 ) 
CRNNBRl 1,9 ) 
CRNNBRl 1,10 ) 
CRNNBRl 1,11 ) 
CRNNBRl 1,12 ) 



IMAX + PLNMAX 
1 
1 

I + 1JMAX-1)*IMAX ♦ PLNMAX 
1 
1 

I+l + PLNMAX 
1 
1 

= I ♦ IMAX ♦ PLNMAX 
= 1 
= 1 



CRmBRl L,l) 
CRNNBRl L,2 ) 
CRmBRl L, 3 ) 



I-l ♦ PLNMAX 

L 

L 






02830 






CRNf'BR(L,^) = I + ( JMAX-1>»IMAX 


♦ PLNMAX 


02840 






CRNNBR(L,5) = L 




02850 






CRNNBR(L,6) = L 




02860 






CRNNBR(L,7) = 1 ♦ PLNMAX 




02870 






CRNfJBRt L ,8 ) = L 




02880 






CRNTJBR( L,9) = L 




02890 






CRNNBR(L,10) = I + IMAX ♦ PLNMAX 


02900 






CRmBR(L>ll) = L 




02910 






CRNNBR(L>12) = L 




02920 


I 


= I MAX 




02922 


J 


= UMAX 




02924 


K 


= 1 






02930 






L = PLNMAX 




02940 






CRNNBR(L,1) = I-l ♦ (J-1)«MAX 


+ PLNMAX 


02950 






CRNNBR( L,2) = L 




02960 






CRNNBR( L,3 ) = L 




02970 






CRNraR(L,4) = I + (J-2)i^IMAX + 


PLNMAX 


02980 






CRNT4BR(L,5) = L 




02990 






CRNNBR( L,6 ) = L 




03000 






CRNNBR(L,7) = 1 +U-1 )i^IMAX + 


PLNMAX 


03010 






CRNr4BR(l,8) = L 




03020 






CRNNBR(L,9) = L 




03030 






CRNNBR(L,10) = I ♦ PLNMAX 




03040 






CRNNBR(L,11) = L 




03050 






CRNNBR(L,12) = L 




03060 


I 


= 1 






03070 




L = 


1 ♦ ( J-1 )^^IMAX 




03080 






CRNNBR(L,1) = IMAX + (J-D^^IMAX 


+ PLNMAX 


03090 






CRNNBR( L,2) = L 




03100 






CRNNBR(.L,3) = L 




03110 






CRNNBR(L,4) = I + (J-2)*IMAX + 


PLNMAX 


03120 






CRNNBR(L,5) = L 




03130 






CRNNBR(L,6) = L 




03140 






CRNNBR( L,7)=I+1 ♦ (J-1)*IMAX + 


PLNMAX 


03150 






CRNraR(L,8) = L 




03160 






CRNNBR( L,9) = L 




03170 






CRNNBR(L,10) = I + PLNMAX 




03180 






CRNNBR(L,11) = L 




03190 • 






CRNNBR( L,12 ) = L 





c 

c 

c 

c 

c 

c 



THE FOLLOWING CODE IS USED TO 
AN OUTPUT FILE. 



SEND THE DEVELOPED ARRAYS TO 



03193 

03195 



0319^ 

03196 

03200 

03210 

03220 

03230 

03240 

03250 

03260 

03270 

03280 



03282 



PRINT 3193, 'THIS IS THE LATICE SITES VS PLANES' 
FORMAT I '1' ,47X,A34) 

PRINT 3195, 'SITE t' , 'PLANE' 

FORMAT ( ,10X,A8,3X,A5) 

DO 3194 L = 1,LMAX 
PRINT 3196, L, 2PLN( U 
COtiTINUE 

FORMAT (' ' ,10X,I8,5X,I3 ) 

PRINT 3210, ' NEAREST NEIGHBOR TABLE ' 

FORMAT ( •!' ,53X,A24) 

PRINT 3230, 'L' ,1,2,3,4,5,6,7,8,9,10,11,12 

FORMAT ( '-' ,4X,A1,9X,I1,7X,I1,7X,I1,7X,I1,7X,I1,7X,I1,7X,I1, 
7X,I1,7X,I1,6X,I2,6X,I2,6X,I2) 

LMAXl = LMAX - PLNMAX 
DO 3270 L = 1, LMAXl 

PRINT 3280, L,(CRNNBR( L,R), R=l,12) 

CONTINUE 

FORMAT ( ' ' ,1318) 

DO 3287 K = 1,KMAX 

PRINT 3284, 'THIS IS PLANE NUMBER', K 
DO 3282 J = 1,JMAX 

PRINT 3286, ( CRYST3( I , J ,K ) , I = 1,IMAX) 
CONTINUE 
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0328 ^ 

03286 

03287 

03288 

03289 



03290 

03291 



03292 

03300 

03301 



03310 

03311 
03320 



FORMAT ( *1' ,5X,A20,3X,I3) 

FORMAT ( * * ,414) 

CONTINUE 

WRITE ( 1,3289 )IMAX ,IMAX1 , JMAX , JMAXl ,KMAX,KMAX1 , PLNMAX , LMAX , LMAXl , 
L ISEED1,ISEED2,MTMSTP,T1,AVGTM,PLSTM,0UTPUT 

FORMAT! • *,918, ‘ *,2116,* *,518) 

DO 3290 L = 1,LMAX 
WRITE 1 1,3291 )CRYST( L) 

CONTINUE 
FORMAT ( * * ,18) 

DO 3300 L = 1, LMAXl 
DO 3292 R = 1,12 

WRITE (1,3301) CRNNBR(L,R) 

CONTINUE 

CONTINUE 

FORMAT ( * * ,318) 

DO 3310 L = 1,LMAX 

WRITE ( 1,3311 )ZPLN(L) 

CONTINUE 

FORMAT ( * * ,18) 

STOP 

END 



APpn::Di:-: c 

'riFPL'SF. ' CODE LISTING 
The following is the c o n p u t e r code for 
s i nu 1 a t i o n n a in e d 'DIFFUSE'. 
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t 'a c diffusion 



•DIFFUSE* 

MARK R. POLNASZEK 
JUNE 1986 



C 
C 

c 
c 
c 
c 

C THE CODE ‘DIFFUSE* HAS BEEN DEVELOPED TO SIMULATE THE DIFFUSION OF 
C AN ATOM INTERSTITIALLY THROUGH A FACE-CENTERED-CUBIC CRYSTAL LATTICE 

C THIS LATTICE IS AN INPUT ARRAY THAT WAS DEVELOPED IN ANOTHER PROGRAM 

C AND IMPUTED ALONG WITH EACH SITES NEAREST-NEIGHBOR, AND EACH SITES 
C PLANAR POSITION. THE METHOD THAT 'DIFFUSE* USES TO MOVE THE 

C DIFFUSING ATOMS THROUGH THE CRYSTAL ARRAY IS BY A MONTE CARLO METHOD 

C THE SIMULATION THEN DEVELOPS THE DIFFUSION COEFFICIENT BY TICK'S 
C FIRST AND SECOND LAW. 

C 

C 

C 

C THE FOLLOWING DEFINES THE VARIABLES AS REAL OR INTEGER AND 
C DIMENSIONS THE VARIABLE ARRAYS* MEMORY ALLOCATION. 

C 

C 

C 

INTEGER PLNMAX, CRYST( 1300 ) ,CRYST3I 08,08 ,16 ) ,CRNIsBR( 1300 ,12 ) ,R , 

1 MAVGTM, ZPLN( 1300),AVGTM,TMSTP2, CONS( 1100 ,16 ) ,T1 ,MTMSTP ,TMSTP , 

2 AVGTMl ,AVGTM2,AT1 ,AVGTM3 ,ATMSTP , OUTPUT 

REAL*4 RDMK 1300 ),RDM2( 1300 ),AC0NS1 210,16), AFLUX( 210 , 16 ) , 

1 ADERCK 210,16) ,ADERC2( 210 ,16 ) ,ADERT( 210 , 16 ) ,DIFC01( 210 ,16 ) , 

2 DIFC02( 210,16), FLUX( 1100 ,16 ), PLNMXl ,CV0L( 1100 ,16 ) , 

3 DERC1(1100,16),DERC2(1100,16),DERT11100,16), ACVOU 210 ,16 ) , 

4 DERDC(210,16),DERT2(210,16) 

READ! 5, ) IMAX,IMAX1 ,JMAX,JMAX1 ,KMAX,KMAX1 ,PLmAX ,LMAX,LMAX1 , 

1 ISEED1,ISEED2,MTMSTP,T1,AVGTM,TMSTP2, OUTPUT 

DO 10 L = 1,LMAX 
C 
C 

C THE FOLLOWING CODE READS IN THE REQUIRED INPUT AND ECHOS 
C THE INPUT INTO AN OUTPUT FILE FOR INSPECTON. 

C 
C 

00010 



00015 

00020 



00030 



00035 

000^0 



000^5 

00050 

00055 

00056 
00060 

1 

00070 



READ (5,*) CRYST(L) 

CONTINUE 

DO 20 L = 1,LMAX1 

DO 15 R = 1,12 

READ i5y^) CRNNBR(L,R) 

CONTINUE 

CONTINUE 

DO 30 L = 1,LMAX 

READ (5,^^) ZPLN(L) 

CONTINUE 

IF (MTMSTP .EQ. 0) GOTO 55 

DO ^0 TMSTP = 1, MTMSTP 

DO 35 K = 1,KMAX1 

READ 15,^) FLUX( TMSTP ,K) 

CONTINUE 

CONTINUE 

DO 50 TMSTP = 1, MTMSTP 

DO ^5 K = 1,KMAX1 

READ iS,^) CONS( TMSTP, K ) 

CONTDnIUE 

CONTIhJUE 

IF (OUTPUT .NE. 1) GOTO 300 

WRITE (6,60) 'IMPUTED VALUES* 

FORMAT ( *1* ,33X,A14) 

WRITE (6,70) ’IMAX','IMAX1*,*JMAX*,'JMAX1*,'KMAX','KMAX1*, 
•PLNMAX* , 'LMAX* , * LMAXl ' , 'TIME AVG * , 'TMSTP 2* 
FORMAT ( •-* ,11A10 ) 

WRITE (6,80) IMAX,IMAX1,JMAX,JMAX1,KMAX,KMAX1, PLNMAX, LMAX, 
LMAXl ,AVGTM , TMSTP 2 
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1 



00080 FORMAT (’ ',11110) 

WRITE (6,90) 'RANDOM (3ENERATOR SEEDS' 

00090 FORMAT ('-' ,29X,A22) 

WRITE (6,0100) ISEED1,ISEED2 
00100 FORMAT ( ' ',10X,2I16) 

WRITE (6,0110) 'TIME CONSTANTS' 

00110 FORMAT ( ' - ' ,33X,A14 ) 

WRITE (6,0120) 'START TIME', 'END TIME' 

00120 FORMAT ('-' ,10X,2A10) 

WRITE (6,0130) T1,MTMSTP 
00130 FORMAT (' ',10X,2I10) 

WRITE (6,0140) 'INPUTED CRYSTAL* 

00140 FORMAT ('!' ,32X,A15) 

WRITE (6,145) 'LATICE ATOM - 2 ',' DIFFUSING ATOM - 1', 
1 'VACANCY - 0* 

00145 FORMAT ('-',3A25) 

WRITE (6,146) 'SITE', 'ATOM TYPE* 

00146 FORMAT ('-' ,10X,2A10 ) 

DO 150 L = 1,LMAX 

WRITE (6,160) CRYST( L ) 

00150 CONTINUE 

00160 FORMAT (' ',10X,2I10) 

WRITE (6,170) 'NEAREST NEIGHBOR TABLE’ 

00170 FORMAT ('!' ,29X,A22 ) 

WRITE (6,0175) ' L ' ,1 ,2 ,3 ,4,5 ,6 , 7 ,8,9,10 ,11 ,12 
00175 FORMAT ( ' - ' , A8 ,1218 ) 

DO 180 L = 1,LMAX1 

WRITE (6,190) L,(CRNNBR( L,R), R=l,12) 

00180 CONTINUE 

00190 FORMAT (' ',1318) 

C 

C 

C THE CODE BELOW SETS UP CONSTANTS THAT ARE USED REPEATEDLY 
C THROUGHOUT THE REMAINDER OF THE PROGRAM. 

C 

C 

00300 PLNMXl = 2/FLOAT( PLNMAX ) 

T1 = T1 + MTMSTP 

MTMSTP = MTMSTP + TMSTP2 
AVGTMl-= AVGTM - 99 
ATI = INT(FL0AT(T1)/AVGTM)+1 
MAVGTM = INT( FLOAT( MTMSTP )/AVGTM ) 

C 

C 

C THE BELOW CODE IS THE BEGINNING OF THE MAJOR DO LOOP. THIS 

C DO LOOP COMPLETES ALL THE REQUIRED TIMESTEPS FOR THE PROGRAM. 

C ALL CALCULATIONS OF THE SIMULATION ARE COMPLETED WITHIN THIS 
C DO LOOP . 

C 

C 

DO 1050 ATMSTP = ATI , MAVGTM 
C 
C 

C THE FOLLOWING DO LOOP COUNTS ALL DIFFUSING ATOMS IN EACH 
C PLANE. THIS IS THE CONCENTRATION OF EACH PLANE. THE DO LOOP 

C ALSO CALCULATES THE CONCENTRATION PER UNIT VOLUME OF EACH 

C CRYSTAL PLANE . 

C 

C 

DO 1001 TMSTP = 1, AVGTM 

DO 301 L = 1,LMAX 
K = ZPLN( L ) 

IF (CRYST(L) .EQ. 1) THEN 
C0N3( TMSTP ,K ) =CONS( TMSTP ,K ) ♦ 1 
CVOL( TMSTP, K) = CVOL( TMSTP ,K ) ♦ PLNMXl 
END IF 

00301 CONTINUE 
C 



c 

C THE FOLLOWING CONSTANT IS DETERMINED FOR USE IN THE OUTPUT FILE. 

C 

C 

AVGTM2 = ( ATMSTP^^AVGTM)-( 100 -TMSTP ) 

C 

C 

C THE FOLLOWING IF STATEMENT SHORTENS THE NUMBER OF SITES THAT THE 

C SIMULATION WILL LOOK AT IF THE STATEMENT IS TRUE, AND ASSIGNS A 

C CONSTANT NUMBER OF SITES IF THE STATEMENT IS FALSE. 

C 

C 

IF (AVGTM2.LE. KMAXl ) THEN 
INMBR = AVGTM2^ PLNMAX 
GOTO 310 
END IF 

INMBR = LMAXl 
C 
C 

C THE FOLLOWING SUB-ROUTINE CALLED FROM THE IMSL LIBRARY SETS UP 
C TWO ARRAYS OF UNIFORM RANDOM NUMBERS. 

C 

C 

00310 CALL SRND( ISEEDl ,RDM1 , INMBR ,2 ,0 ) 

00311 CALL SRND(ISEED2,RDM2, INMBR, 2,0) 

C 

C 

C THE DO 999 DO LOOP LOOKS AT EACH SITE IN THE CRYSTAL ARRAY 
C TO DETERMINE IF A DIFFUSING ATOM IS IN THAT SITE. IF THERE 

C IS A DIFFUSING ATOM IN THAT SITE THE CODE MOVES IT TO NEW 

C SITE IF IT CAN. 

C 

C 

00312 DO 999 LRDM = 1, INMBR 

00313 L = INT(INMBR*RDM1(LRDM)) ♦ 1 

00314 R = INT (12 * RDM2(LRDM)) ♦ 1' 

LI = CRNNBR(L,R) 

C 

C 

C THE FOLLOWING IF STATEMENT DETERMINES IF THE SITE IS OCCUPIED 
C BY A DIFFUSING ATOM. IF THE STATEMENT IS FALSE THE CODE ALLOWS 
C THE DO LOOP TO COrTTINUE. 

C 

C 

IF (CRYST(L) .NE. 1) GOTO 999 
C 
C 

C THE FOLLOWING DO LOOP IS USED TO SEND FIRST PLANE ATOMS TO 
C THEIR SPECIAL CODE. 

C 
C 

00316 
C 
C 

C THE 
C 
C 

00317 

1 

00318 
C 
C 

C THE 
C 
C 

00320 



IF (2PLN(L) .EQ. 1) GOTO 320 



NEXT DO LOOP MOVES THE ATOM TO ITS NEW SITE . 



IF (CRYST(Ll) .EQ. 0) GOTO (400,400,400,400,500,500,500, 

500,600,600,600,600) R 

GOTO 999 



FOLLOWING IF STATEMENT MOVES THE FIRST PLANE ATOMS. 



IF (CRYST(Ll) .EQ. 0) THEN 
FLUX(TMSTP,1) =FLUX(TMSTP,1) + PLNMXl 
CRYST(Ll) = CRYST(L) 



END IF 
GOTO 999 
C 
C 

C THE FOLLOWING CODE MOVES ALL ATOMS OTHER THAN FIRST PLANE 
C IF THE NEAREST-NEIGHBOR SELECTED WAS IN THE FORWARD PLANE. 
C 
C 



00400 


K = ZPLN(L) 




00401 


FLUX( TMSTP, K) = FLUXl TMSTP ,K ) 


4 PLNMXl 


00403 


IF (ZPLN(Ll) .EQ. KMAX) GOTO 


410 


00404 


CRYST(Ll) = CRYST(L) 




00410 


CRYST(L) = 0 






GOTO 999 





c 

c 

C THE FOLLOWING CODE MOVES the aiom if if is fo sfay in fhe same 
C PLANE. 

C 

C 

00500 CRYST(Ll) = CRYST( L) 

CRYST(L) = 0 
GOTO 999 
C 
C 

C THE FOLLOWING CODE IS USED IF THE ATOM IS TO MOVE INTO THE 

C PLANE BEHIND THE DIFFUSING ATOMS INITIAL POSITION. 

C 

C 

00600 K = ZPLN(Ll) 

FLUX(TMSTP,K) = FLUX! TMSTP >K ) - PLNMXl 
CRYST(Ll) = CRYST(L) 

CRYST(L) = 0 
00999 CONTINUE 

01001 COrTTINUE 
C 

C 

C THE FOLLOWING DO LOOP FINDS THE CONCENTRATION FOR THE END OF THE 
C LAST TIMESTEP. 

C 

C 

DO 1002 L = 1,LMAX 
K = 2PLN(L) 

IF (CRYST(L) .EQ. 1) THEN 
CONS( AVGTM+1,K) = CONS( AVGTM+1 ,K ) ♦ 1 
CVOU AVGTM+1,K ) = CVOU AVGTM+1 ,K ) + PLNMXl 
END IF 

01002 CONTINUE 
C 

C 

C THE FOLLOWING DO LOOPS FIGURE THE DERIVATIVES NEEDED TO 
C DETERMINE THE DIFFUSION COEFFICIENTS. 

C 

C 

DO 1004 TMSTP = 1,AVGTM 
DO 1003 K = 2,KMAX1 

DERCK TMSTP, K ) =( CVOU TMSTP ,K+1 ) - CVOL( TMSTP ,K-1 ) )/2 

01003 CO^>mNUE 

01004 CONTINUE 

DO 1006 TMSTP = 1,AVGTM 
DO 1005 K = 3,KMAX1-1 

DERC2( TMSTP ,K ) = ( DERCK TMSTP ,K + 1 )-DERCl( TMSTP ,K-1 ) )/2 

01005 CONTINUE 

01006 CONTINUE 

DO 1008 TMSTP = 1,AVGTM 
DO 1007 K = 1,KMAX1 

DERT( TMSTP ,K ) = CVOL( TMSTP ♦! ,K ) - CVOU TMSTP, K) 

01007 CONTINUE 



7o 



01008 

C 

c 
c 
c 
c 
c 



COfTTINUE 



THE FOLLOWING DO LOOP AVERAGES THE DETERMINED VALUES AS NEEDED 
TO DETERMINE THE DIFFUSION COEFFICIENT, 



01010 

01020 

C 

C 

C 

C 

C 

c 



DO 1020 TMSTP = 
DO 1010 K 
ACONSI ATMSTP,K) 
AFLUXJ ATMSTP,K) 
ADERCH ATMSTP^K) 
ADERC2< ATMSTP,K ) 
ADERT( ATMSTP,K) 
ACVOLl ATMSTP,K) 
CONTINUE 
CONTINUE 



1,AVGT>1 
= 1,KMAX1 
ACONSl ATMSTP,K) + 
AFLUXf ATMSTP,K) ♦ 
= ADERCK ATMSTP,K) 
= ADERC21 ATMSTP,K) 
ADERTl ATMSTP,K) ♦ 
ACVOL(ATMSTP,K) ♦ 



FLOAT( CONSl TMSTP ,K ) )/AVGTM 
F LUX 1 TMSTP, K) /AVGTM 

♦ D E RC 1 1 TMSTP , K ) /AVGTM 

♦ DERC2 1 TMSTP ,K ) /AVGTM 
DERT( TMSTP, K) /AVGTM 

CVOL( TMSTP, K )/AVGTM 



THE FOLLOWING DO LOOP DETERMINES THE DIFFUSION COEFFICIENT 
BY BOTH THE FIRST LAW ( DKOl ) AND SECOND LAW (DIFC02) METHOD. 



DO 1030 K = 3,KMAX1-1 

DIFCOH ATMSTP,K) = -1^^ AFLUX( ATMSTP ,K )/ADERCl( ATMSTP ,K ) 
DIFC02< ATMSTP,K) = ADERTf ATMSTP ,K )/ADERC21 ATMSTP ,K ) 
01030 CONTINUE 

C 
C 

C THE FOLLOWING DO LOOP DETERMINES THE FIRST DERIVATIVE OF THE 
C DIFFUSION COEFFICIENT AND BY FICK*S SECOND LAW f NON-LINEAR) 

C THE TIME DERIVATIVE OF THE CONCENTRATION. 

C 

C 



DO 1031 K = 3,KMAXl-2 

DERDCf ATMSTP,K) = DIFCOK ATMSTP ,K + 1 ) - DIFCOll ATMSTP ,K-1 ) 
DERT2(ATMSTP,K) = DERDC( ATMSTP ,K j^ADERCK ATMSTP ,K ) ♦ 

1 DIFCOK ATMSTP ,K )*ADERC2f ATMSTP ,K ) 

01031 CONTINUE 
C 

C 

C THE FOLLOWING DO LOOP CLEARS THE CONTENTS OF THE NAMED ARRAYS. 

C 

C 

IF (ATMSTP .EQ. MAVGTM ) GOTO 1050 
DO 1033 TMSTP = 1,AVGTM+1 
DO 1032 K = 1,KMAX 

CONS( TMSTP, K ) = 0 
FLUX( TMSTP, K ) = 0 
DERCK TMSTP, K ) = 0 
DERC2( TMSTP, K) = 0 
DERT( TMSTP, K ) = 0 

CVOL( TMSTP, K ) = 0 

01032 CONTI rAJE 

01033 CONTINUE 

01050 CONTINUE 
C 

C 

C THE FOLLOWING CODE WRITES THE GENERATED VALUES OF THE SIMULATION 
C TO AN OUTPUT FILE FOR ANALYSIS. 

C 

C 

GOTO ( 1051, 1191, 1241), OUTPUT 

01051 WRITE (6,1060) *OUTPUTED VALUES* 

01060 FORMAT ( * 1 * , 33X , A15 ) 

WRITE( 6,1070) 'IMAX* , *IMAX1* , ' UMAX* , *JMAX1* , *KMAX* , *KMAX1* , 
1 'PLNMAX* , 'LMAX* , * LMAXl * 

01070 FORMAT ( * - ' , 9A8 ) ' 



1 

01080 

01090 

01100 

OHIO 

01120 

01130 

011^0 

1 

Oll^S 

011^6 



OllSO 

01160 

01170 

0117S 



01180 

01190 

01191 



1 

2 

01200 

01210 

01220 

01230 

01240 

01241 



1 

2 



1 

2 

3 

01340 

01350 

01360 

01370 

01380 



I 

WRITE C 6 ,1080 )IMAX , IMAXl , JMAX , JMAXl ,KMAX ,KMAX1 , PLNMAX , LMAX , I 

LMAXl 

FORMAT (* *,918) , 

WRITE( 6,1090) 'RANDOM GENERATOR SEEDS* I 

FORMAT ( *-* ,29X,A22) . 

WRITE (6,1100) ISEED1,ISEED2 
FORMAT (• *,10X,2I16) 

WRITE (6,1110) 'TIME CONSTANTS* 

FORMAT ( *-* ,33X,A14) | 

WRITE (6,1120) 'START TIME*, *EN0 TIME* j 

FORMAT ( *-* ,10X,2A10 ) 

WRITE (6,1130) T1,MTMSTP 
FORMAT (* *,10X,2I10) 

WRITE (6,1140) 'IMPUTED CRYSTAL* 

FORMAT ( *1* ,32X,A1S) 

WRITE (6,1145) 'LATICE ATOM - 2 *,* DIFFUSING ATOM - 1*, 

'VACANCY - 0* 

FORMAT ( *-* ,3A2S ) 

WRITE(6,1146) 'SITE*, 'ATOM TYPE* 

FORMAT ( *-' ,10Xi2A10) 

DO 1150 L = 1,LMAX 

WRITE( 6,1160 ) L, CRYST( L ) 

CONTINUE 

FORMAT ( * * ,10X,2I10) 

WRITE( 6,11,70 ) 'NEAREST NEIGHBOR TABLE* 

FORMAT ( *1* ,29X,A22) 

WRITE (6,1175) *L* ,1,2,3,4,5,6,7,8,9,10,11,12 
FORMAT ( *-* ,A8,12I8) 

DO 1180 L = 1, LMAXl 

WRITE(6,1190) L,(CRNNBR( L,R), R=l,12) 

CONTINUE 

FORMAT (* *,1318) 

DO 1210 TMSTP = AVGTM1,AVGTM 

AVGTM3 = ( MAVGTM^AVGTM ) - ( 100 -TMSTP ) 

WRITE (6,1220) 'PARAMETERS FOR TIMESTEP t*, AVGTM3 
WRITE (6,1230) 'PLANE*, *CON T- *,' C/VOL T- ',' FLUX * , 

•1ST 0ERV*,'2ND DERV*,*TIME DERV, *CON T+*, 'C/VOL T+ * 

DO 1200 K = 1,KMAX1 

WRITE (6,1240) K ,CONS( TMSTP ,K ) ,CVOL( TMSTP ,K ) , 

FLUX( TMSTP ,K ) , DERCK TMSTP ,K ) ,DERC2( TMSTP ,K ) ,DERT( TMSTP ,K ) , 

CONS(TMSTP+l,K), CVOL( TMSTP+1 ,K ) 

CONTINUE 

CONTINUE 

FORMAT ( *1* ,A23,I6) 

FORMAT ( *-* ,9A10 ) 

FORMAT (• ' ,2110, SFIO. 5,110, FIO. 5) 

DO 1350 ATMSTP = 1,MAVGTM 

WRITE( 6,1360 ) 'AVERAGE PARAMETERS FOR TIMESTEPS*, 

( ( ATMSTP-1 )*AVGTM ) + l, *TO' , ATMSTP^AVGTM 
WRITE (6,1370) * PLANE ','AVG CONC','AVG C/VOL ',*AVG FLUX*, 

*1ST DER*, *2ND DER*, *TIME DER*,*DIF COEF 1*, 

'DIF COEF 2*, *DIF CO DER*, *TIME DER ?* 

DO 1340 K = 1,KMAX1 

WRITE (6,1380 ) K ,ACONS( ATMSTP ,K ) ,ACVOL( ATMSTP ,K ) , 

AFLUX ( ATMSTP ,K ),ADERC1( ATMSTP ,K ) ,ADERC2( ATMSTP,K ), 

ADERT ( ATMSTP ,K ) ,0IFCOl ( ATMSTP ,K ) ,0IFCO2( ATMSTP ,K ) , 

DERDC( ATMSTP, K ),0ERT2( ATMSTP, K ) 

CONTINUE 

CONTIF4UE 

FORMAT ( '1' ,A31,IS,A2,IS ) 

FORMAT ( *-' ,11A11 ) j 

FORMAT (* ' ,I11,10F11.S) I 

STOP , I 

END 



I 







:i:ix B 


SA.hPLE 


lEPUT 


FOR 'BIESr.T' 



The following is the itens that need to be inputed i 
the code ' B I F S E T ' . The n e t h o d to input this data is to 
it in the sane order as shov/n without the titles. There 
no for r. at to tlie length of each data, but insure that tli 
ere two b 1 a n !■: s b e t v; e e n each input. 



1 1 o 
; u t 
i s 
: r e 



79 



IMAX 


JMAX 


KMAX 


MTMSTP 


T1 


ISEEDl 


ISEED2 


04 

AVGTn 

100 


04 

OUTPUT 

1 


06 


00 


01 


447586930 


88475732 



t 




f r o ir. 
This 
T e t 
i n I) u t 



sa;;pli: fcp 'pirF.us;-. 

he f o 1 1 o i n 2 input for the code ' ?• 1 1’ r U S * i s c r o a t c- (' 
the output file, FILZ FTOlFOOl, fror. t:‘.e code 'jTFSr.T'. 
is accomplished by renaming the file as I b’ ? 0 L 2 ''AT.!, 
itles seen in the sample are not a part of the actual 



r; i 



Sample input of the constants required for ‘diffuse* , 

The sample input contains titles that would not be in the 
actual input data* 



IMAX 


IMAXl 


JMAX 


JMAXl 


KMAX 


KMAXl 


PLNMAX 


LMAX 


LMAXl 




3 


4 


3 


6 


S 


16 


96 


80 



RANDOM NUMBER GENERATOR SEEDS 
^^7586950 
e8<^75732 



MTMSTP 

0 



T1 

1 



AVGTM 

100 



TMSTP2 OUTPUT 
100 1 



Sample input for the CRYST(L) array 



Z 

1 

Z 

1 

1 

Z 

1 

Z 

z 

1 

z 

0 

z 

0 

0 

z 

0 

z 

z 

0 



z 

0 

z 

0 

z 

z 

0 

z 

0 

0 

z 

0 

z 

z 

0 

z 

0 



'1 n 



20 

1 

1 

29 

1 

1 

18 

1 

1 

21 

1 

1 

17 

2 

2 

30 

2 

2 

• 

m 

63 

60 

61 

52 



input for the CRNNBR(L,R 



) array. 



9 /. 



Sanple inpui for fbe ZPLNCL) array. 



1 

1 



1 

2 

2 



2 



6 

6 



6 



w J 



sa;:ple outpi 

The s a n'. p 1 G output s h o n 
the output u s e ci to c h e c !■: the 
all of tine output generated, 
generated . 



•'■.’DIP r 

JT OF 'DIFSET’ 



in this Appendix is 


the 


s 0 r;. 


G 0 ^ 


code for accuracy. 


T Pi i 


s i s 


Pi 0 t 


but a representation 


0 L 


\j h a 


t is 



■ 



TE « 

1 

2 

S 

4 

5 

i 

7 

8 

f 

10 

11 

12 

15 

K 

15 

U 

17 

18 

1^ 

20 

21 

22 

25 

24 

25 

2i 

27 

28 

2^ 

50 

51 

52 

55 

54 

55 

54 

57 

58 

55 

40 

41 

42 

45 

44 

45 

44 

47 

48 

45 

50 

51 

52 

55 

54 

55 

54 

57 

58 



THIS IS THE LATICE SITES VS PLANES 



PLANE 

1 

I 

I 

1 

I 

1 

I 

1 

1 

I 

1 

1 
1 
1 
1 
1 

2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
2 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 
4 
4 
4 
4 
4 
4 
4 



4 

4 

4 

4 

4 

4 

4 



^7 



nearest neighbor table 





1 


2 


5 


4 


5 


1 


20 


1 


1 


29 


1 


2 


17 


2 


2 


50 


2 


5 


18 


3 


3 


51 


3 


4 


19 


4 


4 


32 


4 


5 


24 


5 


5 


17 


5 


4 


21 


4 


4 


18 


4 


7 


22 


7 


7 


19 


7 


8 


23 


8 


8 


20 


8 


9 


28 


9 


9 


21 


9 


10 


25 


10 


10 


22 


10 


11 


24 


11 


11 


25 


11 


12 


27 


12 


12 


24 


12 


15 


32 


IS 


15 


25 


IS 


14 


29 


14 


14 ' 


24 


14 


15 


30 


15 


15 


27 


15 


14 


51 


14 


14 


28 


14 


17 


54 


45 


54 


57 


32 


18 


53 


44 


55 


58 


29 


19 


34 


47 


54 


59 


50 


20 


55 


48 


53 


40 


31 


21 


40 


33 


58 


41 


20 


22 


37 


54 


59 


42 


17 


25 


58 


55 


40 


45 


18 


24 


59 


54 


57 


44 


19 


25 


44 


57 


42 


45 


24 


24 


41 


58 


43 


44 


21 


27 


42 


59 


44 


47 


22 


28 


45 


40 


41 


48 


23 


29 


48 


41 


44 


55 


28 


50 


45 


42 


47 


54 


25 


51 


44 


45 


48 


55 


24 


52 


47 


44 


45 


34 


27 


55 


52 


41 


50 


55 


48 


54 


49 


42 


51 


54 


45 


55 


50 


45 


52 


55 


44 


54 


51 


44 


49 


54 


47 


57 


54 


49 


54 


57 


54 


58 


55 


50 


55 


58 


55 


59 


54 


51 


54 


59 


54 


40 


55 


52 


55 


40 


55 


41 


40 


55 


58 


41 


40 


42 


57 


54 


59 


42 


57 


45 


58 


55 


40 


45 


58 


44 


59 


54 


57 


44 


59 


45 


44 


57 


42 


49 


44 


44 


41 


58 


43 


50 


41 


47 


42 


59 


44 


51 


42 


48 


45 


40 


41 


52 


45 


49 


48 


77 


44 


49 


44 


50 


45 


78 


47 


70 


41 


51 


44 


79 


48 


71 


42 


52 


47 


80 


45 


72 


45 


55 


72 


45 


70 


75 


52 


54 


49 


44 


71 


74 


49 


55 


70 


47 


72 


75 


50 


54 


71 


48 


49 


74 


51 


57 


74 


49 


74 


77 


54 


58 


73 


70 


75 


78 


53 


59 


74 


71 


74 


79 


54 


40 


75 


72 


73 


80 


55 


41 


80 


75 


78 


45 


40 


42 


77 


74 


79 


44 


57 


45 


78 


75 


80 


47 


58 


44 


79 


74 


77 


48 


59 



7 


8 


9 


10 


11 


12 


18 


1 


1 


21 


1 


1 


19 


2 


2 


22 


2 


2 


20 


3 


3 


23 


3 


5 


17 


4 


4 


24 


4 


4 


22 


5 


5 


25 


5 


5 


23 


4 


4 


24 


4 


4 


24 


7 


7 


27 


7 


7 


21 


8 


8 


28 


8 


0 


24 


9 


9 


29 


9 


9 


27 


10 


10 


30 


10 


10 


20 


11 


11 


31 


11 


11 


25 


12 


12 


32 


12 


12 


50 


IS 


13 


17 


IS 


13 


51 


14 


14 


18 


14 


14 


52 


15 


15 


19 


15 


15 


29 


14 


14 


20 


14 


14 


22 


24 


4 


13 


2 


5 


25 


21 


1 


14 


3 


4 


24 


22 


2 


15 


4 


7 


21 


23 


5 


14 


1 


8 


24 


28 


8 


1 


4 


9 


27 


25 


5 


2 


7 


10 


28 


24 


4 


3 


8 


11 


25 


27 


7 


4 


5 


12 


50 


32 


12 


5 


10 


IS 


51 


29 


. 9 


4 


11 


14 


52 


50 


10 


7 


12 


15 


29 


31 


11 


6 


9 


14 


18 


20 


14 


9 


14 


1 


19 


17 


13 


10 


15 


2 


20 


18 


14 


11 


14 


5 


17 


19 


15 


12 


IS 


6 


58 


40 


20 


29 


18 


21 


59 


57 


17 


50 


19 


22 


40 


38 


18 


51 


20 


25 


57 


59 


19 


52 


17 


24 


42 


44 


24 


17 


22 


25 


43 


41 


21 


18 


25 


24 


44 


42 


22 


19 


24 


27 


41 


45 


23 


20 


21 


28 


44 


48 


28 


21 


24 


29 


47 


45 


25 


22 


27 


30 


48 


44 


24 


25 


28 


51 


45 


47 


27 


24 


25 


52 


34 


54 


32 


25 


50 


17 


55 


33 


29 


24 


51 


18 


54 


34 


30 


27 


52 


19 


55 


35 


51 


28 


29 


20 


54 


54 


54 


45 


54 


57 


55 


55 


55 


44 


35 


38 


54 


54 


54 


47 


54 


59 


55 


55 


55 


48 


55 


40 


58 


40 


40 


53 


58 


41 


59 


57 


57 


54 


59 


42 


40 


58 


58 


55 


40 


45 


57 


59 


59 


54 


57 


44 


42 


44 


44 


57 


42 


45 


43 


41 


41 


58 


43 


44 


44 


42 


42 


39 


44 


47 


41 


45 


43 


40 


41 


48 


50 


52 


48 


41 


44 


33 


51 


49 


45 


42 


47 


34 


52 


50 


44 


43 


48 


55 


49 


51 


47 


44 


45 


54 



4 

1 

2 

5 

< 

s 

i 

7 

8 

10 

11 

12 

15 

1 < 

IS 

14 

50 

51 

52 

2 ^ 

18 

20 

17 

22 

25 

24 

21 

24 

27 

28 

25 

44 

47 

48 

45 

54 

55 

54 

35 

58 

59 

40 

57 

42 

45 

44 

41 

42 

45 

44 

41 

50 

51 

52 

4f 

54 

55 

54 

55 

58 

59 

40 

57 



THIS IS PLANC NUMBER 



I 



2 12 1 
12 12 
2 12 1 
12 12 






THIS IS PLANE NUHBER 



5 



2 0 2 0 
0 2 0 2 
2 0 2 0 
0 2 0 2 
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AP?i::: 


FIX G 






SAMPLE OUTPUT 


or 'DIFFUSE' 






E' ii e f 0 1 1 o i n ^ output is a 


representation of t;ie 


out 


pu 


that 'DIFFUSE' generates to 


use in t li e analysis 


0 f 


th 



diffusion prcoss. 
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OUTPUTED VALUES 



IMAX 


IMAXI 


JMAX 


UMAX! 


KHAX 


KMAXt 


PLNMAX 


LMAX 


lmaxi 


A 


3 


4 


3 


4 


5 


14 


44 


80 



RANDOM GENERATOR SEEDS 
8S426SA2S 18Ub88SU 



TIME CONSTANTS 



START TIME END TIME 
t 100 



0 9 



IMPUTED CRYSTAL 



LATICE atom - 2 D1FFU31MG ATOM - I 



VACANCY - 0 



SITE 

1 

2 

S 

4 

5 
ft 
7 
S 
9 

10 

11 

12 

IS 

1 < 

IS 

1ft 

17 

18 
\9 
20 
21 
22 
25 

24 

25 

26 

27 

28 
29 

50 

51 

52 
55 

54 

55 
5ft 
57 

38 

39 

40 

41 

42 
45 

44 

45 
4ft 

47 

48 

49 

50 

51 

52 
55 

54 

55 
5ft 
57 
SB 

59 

60 
61 
62 
65 
64 



ATOM TYPE 
2 
1 
2 
1 
1 
2 
1 
2 
2 
1 
2 
1 
1 
2 
1 
2 
0 
2 
0 
2 
2 
0 
2 
0 
0 
2 
0 
2 
2 
0 
z 
0 
2 
0 
2 
0 
0 
2 
0 
2 
2 
0 
2 
0 
0 
*> 

0 

2 

0 

2 

0 

2 

2 

0 

2 

0 

0 

2 

0 

2 

2 

0 

2 

0 



Pi P 
- u 



I 

20 

17 

18 

If 

24 

21 

22 

25 

28 

2S 

24 

27 

52 

2f 

50 

51 

54 

55 

56 

5S 

60 

57 

58 

5f 

66 

61 

62 

65 

68 

65 

64 

67 

52 

6f 

50 

51 

54 

55 

56 

55 

40 

57 

58 

5f 

46 

41 

62 

45 

48 

45 

44 

47 

72 

49 

70 

71 

74 

75 

76 

75 

80 

77 

78 

79 



12 

1 

5 

6 

5 

4 

7 

8 

9 

10 

11 

12 

15 

16 

15 

14 

S 

4 

7 

8 

9 

10 

11 

12 

15 

16 

15 

14 

1 

2 

5 

6 

22 

25 

26 

25 

24 

27 

28 

29 

30 

51 

52 

17 

18 

19 

20 

57 

58 

39 

60 

61 

62 

65 

66 

65 

64 

67 

68 

55 

56 

35 

54 



nearest neighbor table 



2 


5 


4 


5 


4 


7 


8 


9 


10 


1 1 


1 


1 


29 


1 


1 


18 


1 


1 


21 


1 


2 


2 


50 


2 


2 


19 


2 


2 


22 


2 


5 


5 


51 


5 


5 


20 


5 


5 


25 


5 


4 


4 


52 


4 


4 


17 


4 


4 


26 


4 


5 


5 


17 


5 


5 


22 


5 


5 


25 


5 


4 


4 


18 


4 


4 


25 


4 


4 


24 


4 


7 


7 


19 


7 


7 


24 


7 


7 


27 


7 


8 


8 


20 


8 


8 


21 


8 


8 


28 


8 


9 


9 


21 


9 


9 


24 


9 


9 


29 


9 


10 


10 


22 


10 


10 


27 


10 


10 


50 


10 


11 


11 


25 


1 1 


1 1 


28 


11 


11 


51 


11 


12 


12 


24 


12 


12 


25 


12 


12 


52 


12 


15 


15 


25 


15 


15 


50 


15 


15 


17 


15 


14 


16 


24 


14 


16 


51 


14 


16 


18 


14 


15 


IS 


27 


15 


15 


52 


15 


15 


19 


15 


14 


14 


28 


14 


14 


29 


14 


14 


20 


14 


65 


54 


57 


52 


50 


22 


24 


4 


15 


2 


64 


55 


58 


29 


51 


25 


21 


1 


14 


5 


67 


54 


59 


50 


52 


26 


22 


2 


15 


6 


68 


55 


60 


51 


29 


21 


25 


5 


14 


1 


55 


58 


61 


20 


18 


24 


28 


8 


1 


4 


56 


59 


62 


17 


19 


27 


25 


5 


2 


7 


5S 


60 


65 


18 


20 


28 


24 


4 


5 


8 


54 


57 


66 


19 


17 


25 


27 


7 


6 


5 


57 


62 


65 


26 


22 


50 


52 


12 


5 


10 


58 


65 


64 


21 


25 


51 


29 


9 


4 


1 1 


59 


66 


67 


22 


26 


52 


50 


10 


7 


12 


60 


61 


68 


25 


21 


29 


51 


1 1 


8 


9 


61 


64 


55 


28 


24 


18 


20 


14 


9 


14 


62 


67 


54 


25 


27 


19 


17 


15 


10 


15 


65 


68 


55 


24 


28 


20 


18 


14 


11 


14 


66 


65 


54 


27 


25 


17 


19 


15 


12 


15 


41 


50 


55 


68 


64 


38 


60 


20 


29 


18 


42 


51 


56 


65 


67 


59 


57 


17 


50 


19 


45 


52 


55 


64 


68 


60 


58 


18 


51 


20 


46 


69 


54 


67 


65 


57 


59 


19 


52 


17 


69 


56 


57 


54 


56 


62 


66 


24 


17 


22 


50 


55 


58 


35 


55 


65 


61 


21 


18 


25 


51 


54 


59 


36 


34 


66 


62 


22 


1 9 


26 


52 


55 


40 


55 


55 


61 


65 


25 


20 


21 


55 


58 


41 


60 


58 


64 


68 


28 


21 


24 


56 


59 


42 


57 


59 


47 


65 


25 


22 


27 


55 


40 


45 


58 


60 


68 


64 


24 


25 


28 


54 


57 


46 


59 


57 


45 


67 


27 


26 


25 


57 


42 


69 


66 


62 


56 


54 


52 


25 


50 


58 


45 


50 


61 


65 


55 


55 


29 


24 


51 


59 


46 


51 


42 


66 


54 


56 


50 


27 


52 


40 


41 


52 


65 


61 


55 


55 


31 


28 


29 


77 


44 


49 


46 


42 


56 


54 


54 


65 


34 


78 


47 


70 


41 


45 


55 


55 


55 


64 


35 


79 


48 


71 


42 


46 


54 


56 


56 


67 


34 


80 


45 


72 


45 


41 


55 


55 


55 


68 


33 


45 


70 


75 


52 


50 


58 


40 


60 


55 


38 


44 


71 


76 


69 


51 


59 


57 


57 


54 


59 


47 


72 


75 


50 


52 


40 


58 


58 


55 


60 


48 


49 


74 


51 


69 


57 


59 


59 


54 


37 


49 


76 


77 


54 


54 


42 


44 


66 


57 


62 


70 


75 


78 


55 


55 


45 


41 


61 


58 


65 


71 


74 


79 


54 


54 


44 


42 


62 


59 


64 


72 


75 


80 


55 


55 


41 


45 


65 


60 


61 


75 


78 


45 


40 


58 


50 


52 


68 


61 


64 


76 


79 


44 


57 


59 


51 


49 


65 


62 


47 


75 


80 


47 


58 


40 


52 


50 


64 


65 


68 


74 


77 


48 


59 


57 


49 


51 


67 


64 


65 
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PARAMETERS FOR TIMESTEP 100 



PLANE 


CON T- 


C/VOL T- 


FLUX 


1ST DERV 


2ND DERV 


TIME DERV 


CON T* 


C/VOL T* 


1 


8 


1.00000 


0.00000 


O.OOOOO 


0.00000 


0.00000 


8 


1.00000 


Z 


6 


0.75000 


0.00000 


-0 . 12500 


0.00000 


0.00000 


6 


0.75000 


5 


6 


0.75000 


0.12500 


-0.51250 


-0.12500 


-0.12500 


5 


0.62500 


< 


1 


0.12500 


0.00000 


-0.57500 


0.12500 


0.12500 


2 


0.25000 


5 


0 


0.00000 


0.00000 


-0.06250 


0.00000 


0.00000 


0 


0.00000 
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AVERAGE parameters FOR TIMESTEP ITO 100 



PLANE 


AVG CONC 


AVG C/VOL 


AVG FLUX 


1ST DER 


2ND DER 


time DER 


DIF COEF 1 


DIF COEF 2 DIF CO DER 


TIME DER 7 


1 


7. ^^999 


1.00000 


0 . 08500 


0.00000 


0.00000 


0.00000 


0.00000 


0.00000 0.00000 


0.00000 


2 




0.7S125 


0.07750 


-0.23812 


0.00000 


0.00750 


0.00000 


0.00000 0.00000 


0.00000 


5 


4. 18^97 


0.52375 


0.07125 


-0.2U25 


0.03187 


0.00625 


0.32948 


0.19608 0.39427 


-0.07474 


4 


2.5A990 


0. 3187'.« 


0.06875 


-0.17437 


0.02844 


0.00250 


0.39427 


0.08791 0.00000 


0 . 00000 


S 


1.39999 


0.17500 


0.06875 


-0,15937 


0.00000 


0.00000 


0.00000 


0.00000 0.00000 


0.00000 
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